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Abstract 

The  locations  of  Hopf  bifurcation  points  associated  with  the  viscous,  incom¬ 
pressible  flow  about  a  NACA  0012  airfoil  with  structural  coupling  are  computed  for 
very  low  Reynolds  numbers  (<  2000).  A  semi-implicit,  first-order-accurate  time- 
integration  algorithm  is  employed  to  solve  the  streamfunction-vorticity  form  of  the 
Navier-Stokes  equations.  The  formulation  models  the  inclusion  of  simple  structural 
elements  affixed  to  the  airfoil  and  captures  the  resulting  airfoil  motion.  The  equations 
describing  the  airfoil  motion  are  integrated  in  time  using  a  fourth-order  Runge-Kutta 
algorithm. 

The  dissertation  is  divided  into  two  parts.  In  Part  I,  numerical  experiments 
are  performed  in  the  laminar  regime  to  determine  if  the  structural  model  of  the 
airfoil  has  an  effect  upon  the  location  of  the  Hopf  bifurcation  point  when  compared 
with  the  fixed  airfoil.  Results  are  reported  for  a  variety  of  structural  characteris¬ 
tics,  including  variation  of  torsional  and  linear  spring  constants,  inertial  properties, 
structural  coupling,  and  structural  damping.  The  structure  of  the  solution  space  is 
explored  by  means  of  phase  plots. 

In  Part  II,  the  Baldwin-Lomax  turbulence  model  is  implemented  to  model 
turbulent  flow.  A  numerical  effort  is  made  to  predict  the  onset  of  unsteady  flow  and 
the  results  are  compared  to  theory  and  experiment.  A  comparison  is  undertaken 
with  the  compressible  and  incompressible  form  of  the  Navier-Stokes  equations  to 
assess  the  correspondence  of  the  unsteady  lift  evolved  by  pitching  the  airfoil  at  a 
prescribed  frequency  and  amplitude.  Comparisons  are  made  with  the  lift  predicted 
by  the  potential  flow  method  of  Theodorsen. 
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HOPF  BIFURCATION  IN  VISCOUS,  LOW  SPEED  FLOWS 
ABOUT  AN  AIRFOIL  WITH  STRUCTURAL  COUPLING 


I.  Introduction 


1.1  Overview  of  Part  I 

In  the  past  few  years,  interest  in  low  Reynolds  number,  time-periodic  flows 
haw  been  increasing.  Pulliam  [1,  2]  has  performed  numerical  studies  detailing  such 
flows,  demonstrating  such  phenomena  as  period-doubling  bifurcations  and  chaotic 
attractors.  Pulliam  investigated  periodic  flows  using  time  integration  of  the  Navier- 
Stokes  equations,  and  found  a  sequence  of  period-doubling  bifurcations  leading  to 
chaos.  Jackson  [3]  has  investigated  the  onset  of  vortex  shedding  corresponding  to 
time-periodic  body  forces  for  a  variety  of  body  shapes.  Strganac  and  Mock  [4] 
have  numerically  investigated  subsonic  flutter  for  a  finite  wing  using  a  potential  flow 
model,  demonstrating  the  existence  of  a  “flutter  boundary,”  above  which  limit-cycle 
oscillations  persist.  A  substantial  portion  of  Part  I  w as  reported  earlier  [5]  as  part 
of  an  ongoing  research  effort. 

The  development  of  a  time-periodic  flow,  emanating  from  an  equilibrium  flow, 
is  an  example  of  Hopf  bifurcation  [6,  7].  Physically,  Hopf  bifurcation  represents  an 
exchange  of  stability  manifested,  in  the  fluid  dynamics  context,  as  the  formation  of  a 
time-periodic  vortex  street.  The  physically  relevant  processes  occuring  are  the  con¬ 
vection  and  dissipation  of  vorticity.  If  the  dissipative  mechanism  dominates,  vorticity 
will  be  dissipated  before  a  wake  instability  can  be  excited,  and  thus  a  steady  flow 
results.  Conversely,  if  the  convective  forces  dominate,  unsteady  vortex  shedding  en¬ 
sues.  The  Reynolds  number,  representing  the  ratio  of  convective  to  dissipative  forces, 
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is  therefore  physically  relevant  and  serves  as  the  bifurcation  parameter.  The  proper¬ 
ties  of  the  solution  space,  as  the  bifurcation  parameter  changes,  serve  to  classify  the 
bifurcation.  In  this  study,  a  subcritical  bifurcation  is  demonstrated  numerically.  A 
subcritical  bifurcation  occurs  when,  across  a  given  range  of  the  bifurcation  parame¬ 
ter,  three  solution  branches  to  the  equation  set  coexist:  a  stable  equilibrium  branch, 
a  stable  limit-cycle  branch,  and  an  unstable  limit-cycle  branch  connecting  the  stable 
branches  [6].  This  region  is  referred  to  in  the  literature  as  bistable.  Physically,  a  sys¬ 
tem  exhibiting  a  subcritical  bifurcation  will  be  susceptible  to  perturbations.  Seydel 
[6]  offers  a  more  complete  description  of  this  type  of  bifurcation. 

Beran  [8]  has  employed  time  integration  as  a  means  of  determining  and  veri¬ 
fying  the  location  of  Hopf  points  for  a  fixed  airfoil  in  a  two-dimensional  parametric 
space  defined  by  Reynolds  number  and  angle  of  attack.  This  was  accomplished  in 
conjunction  with  a  direct  method  for  computing  the  Hopf  point  [3,  9].  The  current 
effort  represents  a  modification  to  the  time-integration  algorithm  of  Beran  to  allow 
the  numerical  prediction  of  the  Hopf  bifurcation  point  for  an  airfoil  with  a  linear 
and  torsional  spring  affixed.  This  represents  a  simple  structural  model  and  thus  sug¬ 
gests  a  means  to  examine  “low  Reynolds  number”  flutter.  The  centred  question  to 
be  investigated  is  the  coupling  effect  between  the  developing  aerodynamic  flowfield 
and  the  structural  model  of  the  airfoil,  which  is  now  free  to  vibrate  or  oscillate  in 
response  to  the  aerodynamic  forces.  The  basic  assumptions  associated  with  the  fluid 
are  those  of  am  incompressible,  viscous,  two-dimensional  flow.  The  fluid  equations 
are  presented  in  streamfunction-vorticity  form. 

The  structurad  model  of  the  airfoil  allows  motion  with  two  degrees  of  freedom: 
pitch  and  vertical  displacement.  The  equations  of  motion  for  the  aurfoil  incorporate 
linear  and  torsional  spring  constants,  structural  daunping  in  both  axes,  amd  struc¬ 
tural  coupling  effects.  Numerical  experiments  are  undertaken  to  determine  the  effect 
on  the  critical  Reynolds  number,  Re„it,  when  the  inertial  amd  structural  properties 
of  the  aurfoil  are  varied.  A  limited  grid  refinement  study  is  undertaken  to  evalu- 
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ate  the  accuracy  of  the  computed  Strouhal  number  and  Re^t-  The  variation  of 
Strouhal  number  with  the  computational  domain  size  is  also  explored.  In  addition, 
the  structure  of  the  Hopf  bifurcation  is  investigated,  with  emphasis  on  the  impact  of 
the  airfoil  structural  model.  These  investigations  are  undertaken  in  the  laminar  flow 
regime  in  order  to  explore  the  solution  structure  in  the  absence  of  turbulence  and 
unnecessary  amounts  of  artificial  dissipation.  Low-speed  flutter  and  the  associated 
solution  structure  in  the  laminar  regime  have  not  been  investigated  for  viscous  flows. 

1.1.1  The  Influence  of  Turbulence  The  transition  from  laminar  to  turbulent 
flow  has  a  fundamental  and  significant  effect  upon  the  onset  of  unsteady  motion  in  a 
fluid.  The  physical  progression  of  this  transition  is  shown  in  Figure  1  for  the  case  of 
a  circular  cylinder  [10].  At  very  low  Reynolds  numbers  (Re  <  5),  the  flow  is  steady 
and  no  vortices  are  present.  As  the  Reynolds  number  increases,  the  flow  remains 
steady,  but  a  pair  of  standing  vortices  appear  in  the  near  wake.  A  further  increase  in 
Reynolds  number  induces  unsteady  laminar  flow,  manifested  by  asymmetric  vortex 
shedding.  It  was  this  laminar  regime  that  was  explored  for  the  moving  airfoil  in 
Part  I  to  assess  the  impact  of  structural  model  on  the  flow  structure. 

The  transition  to  turbulence  at  larger  Reynolds  numbers  is  the  hallmark  of 
the  next  regime.  As  the  Reynolds  number  increases  to  approximately  3  x  105,  the 
increased  dissipation  provided  by  turbulence  disrupts  the  unsteady  vortex  street  and 
a  quasi-steady  flow  ensues.  As  the  Reynolds  number  increases  beyond  3.5  x  106,  the 
flow  again  becomes  unsteady  with  the  reestablishment  of  the  vortex  street.  The 
transition  to  turbulent  flow  represents  a  natural  demarcation  between  Parts  I  and  II 
of  this  work  as  applied  to  the  prediction  of  flutter  onset.  The  numerical  investigations 
in  Part  I  are  applied  in  the  flow  regimes  indicated  by  the  first  three  entries  in  Figure  1. 
The  numerical  investigation  of  flutter  onset  in  Part  II  incorporates  a  turbulence 
model,  and  is  applied  in  the  flow  regimes  indicated  by  the  final  two  entries  in  Figure  1. 
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w»  <  »  REGIME  OF  UN3EPARATEO  FLOW 


5  TO  15  <  R«  <  40  A  FIXED  PAIR  OF  FOPPL 
VORTICES  IN  WAKE 


40  <  R«  <  90  AND  90  <  R«  <  ISO 
TWO  REGIMES  IN  WHICH  VORTEX 
STREET  IS  LAMINAR 


ISO  <  R«  <  300  TRANSITION  RANGE  TO  TURBU¬ 
LENCE  IN  VORTEX 

300  <  R.  ^  3  X  I05  VORTEX  STREET  IS  FULLY 
TURBULENT 


3  X  IQ5  g  R»  <  3  5  X  108 

LAMINAR  BOUNOARY  LAYER  HAS  UNOERGONE 
TURBULENT  TRANSITION  AND  WAKE  IS 
NARROWER  AND  DISORGANIZED 


3.5  X  106  <  R« 

RE-ESTABLISHMENT  OF  TURBU¬ 
LENT  VORTEX  STREET  . 


Figure  1.  Regimes  of  fluid  flow  across  a  circular  cylinder  [10] 
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1.2  Overview  of  Part  II 

While  the  exploration  of  the  solution  space  presented  in  Part  I  is  of  value 
from  an  analytic  point  of  view,  the  results  do  not  directly  correspond  with  ex¬ 
periment  (typically  performed  in  the  subsonic  to  transonic  regimes  [11,  12],  with 
Re  >  5  x  10s,  where  turbulence  is  a  factor).  The  predicted  onset  of  unsteady 
motion  (laminar  regime)  occurs  at  a  Reynolds  number  much  lower  than  the  exper¬ 
imental  evidence  (turbulent  regime)  indicates.  Therefore,  an  attempt  is  made  to 
capture  experimental  data  by  implementing  a  turbulence  model  to  account  for  these 
effects  at  larger  Reynolds  numbers.  The  turbulence  model  selected  for  this  effort 
is  that  of  Baldwin  and  Lomax  [13]  as  modified  by  Knight  and  Visbal  [14].  A  brief 
discussion  of  the  implementation  of  the  turbulence  model  and  its  applicability  to 
the  present  calculations  is  provided  in  Chapter  4.  An  attempt  is  made  to  capture 
numerically  the  experimental  data  of  Yang  and  Zhao  [11]  and  the  theoretical  and 
experimental  results  of  Theodorsen  [15]. 

A  similar  computational  effort,  on  a  larger  scale,  was  accomplished  by  Gu- 
ruswamy  [16].  He  employed  the  compressible  Navier-Stokes  equations  with  the 
Baldwin-Lomax  turbulence  model  in  conjunction  with  an  aeroelastic  model  of  a  wing- 
body  configuration  and  successfully  simulated  aeroelastic  oscillations  for  a  transonic 
flow  regime  (M0 0  =  0.975).  Although  Guruswamy  achieved  satisfactory  results  using 
the  Baldwin-Lomax  turbulence  model,  and  demonstrated  no  specific  shortcomings, 
he  did  question  the  accuracy  of  the  turbulence  model  in  predicting  the  onset  of  flow 
unsteadiness  due  to  the  quasi-steady  assumptions  applied  in  its  derivation.  Kousen 
and  Bendiksen  [17]  and  Wu  et  al.  [18]  have  also  applied  the  Navier-Stokes  equa¬ 
tions  to  the  numerical  prediction  of  flutter  onset  in  the  high  subsonic  and  transonic 
regimes.  Similar  efforts  applying  the  Euler  equations  have  been  reported  by  Bendik¬ 
sen  and  Kousen  [19]  and  Robinson  et  al.  [20].  Strganac  et  al.  [21]  investigated  the 
numerical  simulation  of  subsonic  flutter  using  a  vortex-lattice  method.  Considerable 
attention  has  been  devoted  to  the  understanding  and  prediction  of  the  onset  of  aero- 
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dynamic  flutter  [12,  15,  22,  23,  24,  25,  26].  In  addition  to  the  computational  efforts 
cited,  wind  tunnel  [22,  23]  and  flight  testing  [12,  24]  in  conjunction  with  theoretical 
analysis  [15,  25]  have  been  the  usual  mode  of  inquiry.  The  advantages  of  a  com¬ 
putational  approach  include  reduced  cost,  easier  application  (once  a  validated  code 
is  available),  and  an  ability  to  explore  experimental  configurations  or  parametric 
studies. 

The  contribution  attempted  in  Part  II  is  twofold.  Firstly,  numerical  evidence 
is  presented  that  the  bifurcation  to  an  unsteady  state  is  subcritical  in  the  turbulent 
regime.  Secondly,  a  basis  is  sought  establishing  the  validity  of  the  numerical  predic¬ 
tion  of  flutter  onset.  In  order  to  substantiate  the  basis  for  numerically  predicting 
flutter  onset,  a  correlation  is  made  between  the  present  numerical  results  and  those 
predicted  by  the  aerodynamic  transfer  function  of  Theodorsen  [15].  A  NACA  0012 
airfoil  is  oscillated  with  a  prescribed  frequency  and  amplitude;  the  computed  lift 
coefficient  is  compared  with  that  predicted  by  Theodorsen ’s  function.  This  is  ac¬ 
complished  using  a  modified  version  of  the  incompressible  code  employed  in  Part  I, 
and  a  Beam  and  Warming  code  authored  by  Visbal  [27],  This  is  the  focus  of  the 
second  part  of  the  present  research  effort.  The  results  are  presented  in  Chapter  5. 

The  modifications  to  the  incompressible  code  are  presented  in  Chapter  4.  A 
time-dependent  metric  transformation  is  applied  to  account  for  the  airfoil  motion 
rather  than  the  approach  outlined  in  Part  I.  The  reasons  for  the  modifications  and 
the  resultant  changes  to  the  boundary  conditions  are  discussed.  Also  provided  in 
Chapter  4  is  an  outline  of  the  development  of  Theodorsen’s  function  and  the  of 
application  of  the  Beam  and  Warming  algorithm. 

Validation  of  the  incompressible  code  employed  in  this  work  is  presented  in  Ap¬ 
pendix  A.  The  implementation  of  the  integration  procedure  for  solving  the  Navier- 
Stokes  equations  was  validated  by  simulation  of  the  flow  about  a  fixed  circular 
cylinder.  The  circular  cylinder  was  chosen  because  a  wide  range  of  experimental 
data  is  available  [28,  29,  30].  In  contrast,  there  is  a  paucity  of  experimented  and 
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numerical  data  available  for  airfoils  over  the  very  low  Reynolds  number  range  con¬ 
sidered.  Therefore,  comparisons  with  experimental  and  other  numerical  data  are  ac¬ 
complished  for  the  Strouhal  number  and  drag  coefficient  for  the  case  of  the  circular 
cylinder.  Good  agreement  was  achieved  with  the  results  of  previous  investigations. 

The  airfoil  equations  of  motion  are  identical  in  form  (second-order  ordinary 
differential  equations)  and  thus  only  the  vertical  axis  was  examined  for  validation 
purposes.  The  validation  was  accomplished  by  comparing  an  exact  solution  to  the 
ordinary  differential  equation  representing  the  vertical  axis  to  the  numerical  solution 
when  a  prescribed  forcing  function  is  applied.  A  phase  lag  equal  to  the  largest 
timestep  used  in  the  airfoil  computations  was  applied  between  the  forcing  functions 
to  examine  the  effect  on  the  respective  responses.  The  method  of  coupling  the 
aerodynamics  and  the  structural  model  was  demonstrated  to  be  accurate. 
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II.  Analysis:  Part  I 


2. 1  Equations  of  Motion 

The  equations  of  motion  and  the  method  of  solution  are  presented  in  this  sec¬ 
tion.  The  modifications  to  the  streamfunction-vorticity  formulation  of  the  equations 
are  described  for  the  case  of  the  rotating  and  translating  airfoil.  Finally,  the  imple¬ 
mentation  of  the  numerical  algorithm  is  discussed,  including  the  application  of  the 
boundary  conditions. 

Equations  are  presented  throughout  in  nondimensional  form.  The  velocity  scale 
is  the  freestream  velocity,  C/,  the  length  scale  is  the  chord  length,  c,  time  is  nondi- 
mensionalized  by  the  aerodynamic  scale,  cfU .  Force  per  unit  mass  (acceleration)  is 
nondimensionalized  by  U2/c. 


In  modifying  the  equations  of  motion  for  a  two-dimensional,  incompressible 
flow,  the  approach  of  Batchelor  [31]  is  applied,  wherein  the  moving  frame,  i.e.,  the 
frame  fixed  to  the  airfoil,  hereafter  referred  to  as  the  a-frame,  is  treated  as  the  com¬ 
putational  reference  frame.  The  acceleration  of  the  a-frame  relative  to  a  true  inertial 
frame  is  rectified  by  inclusion  of  apparent  body  forces.  These  body  forces  (per  unit 
mass,  nondimensionalized)  are  then  added  directly  to  the  momentum  equation  in 
vector  form.  The  components  of  momentum  must  then  be  examined  individually  to 
insure  they  are  correctly  handled  when  transforming  to  the  streamfunction-vorticity 
form  of  the  equations.  An  alternate  approach  to  modeling  the  moving  airfoil,  involv¬ 
ing  a  time-dependent  metric  transformation,  is  presented  in  Part  II. 


The  acceleration  in  the  inertial  system  is  given  by  Greenwood  [32]: 


(PR  <Pf  -  dr  dQ  ^  ~  ^ 

“  =  W  +  ^  +  2fix^+ ^rx’-  +  nx(nxr-), 


(PR  Du 
dt 2  +  Dt 


dft  , 

+  -?rr  +  217  x  u  +  —  x  r  +  fi  x  (fi  x  r). 

at 


(i) 
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where  R(t)  is  the  position  vector  to  the  origin  of  the  a-frame,  r  represents  the  position 
vector  of  a  fluid  element  in  the  a-frame,  and  H(l)  is  the  angular  velocity  of  the  a- 
frame  with  respect  to  the  fixed  frame  (see  Figure  2).  Here  the  term  2  fl  x  u  accounts 
for  the  Coriolis  effect,  dQ/dt  x  r  is  associated  with  tangential  (Eulerian)  acceleration, 
and  fi  x  (fl  x  r)  is  the  centripetal  acceleration,  while  u  =  uex  4  ve2  represents  the 
velocity  in  the  a-frame.  The  velocity  observed  in  the  inertial  system  is  related  to 
this  velocity  by  the  relation  [32] 

4  u  +  n  x  f,  (2) 

or,  in  components  referenced  to  the  a-frame  system, 

V  =  (J?i  4  u  -  ySl)ei  -I-  (R2  4  v  +  xfl)e2 

=  uei  +  ve2.  (3) 

The  velocity  field  observed  in  the  inertial  frame  expressed  in  components  referenced 
to  that  frame  is  given  by 


V  =  (ucosa  4-  v  sin  a)i  +  (  — usina  +  t;coso)_;  , 


(4) 


where  a  =  a(t)  and  a  =  — fl. 

The  acceleration  in  the  a-frame  may  be  equated  to  the  local  force  per  unit 
mass  acting  on  a  fluid  element  to  give  the  momentum  equation.  In  terms  of 
assuming  an  Eulerian  formulation  [31], 


Du 

Dt 


du  _  <Pr 

=  —  +  u  ■  Vu  =  — — . 
df  dt 2 


(5) 


The  applied  forces  can  be  separated  into  the  pr-  s. 
and  the  apparent  body  forces.  If  the  body  forces  an 


.'p  forces,  the  viscous  forces, 
— * 

•  t  oted  by  /  =  /jei  +  /2e2, 
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then  the  components  of  the  momentum  equation  [31]  may  be  written 


du  du  du 
dv  dv  dv 

at+uai+vai 


dp  1  /  d2u  d2u\ 

lh  +  lu{foi  +  lw)~fu 

dp  1  /  d2v  d2t;\ 

dy  +  ~Re  V^2  +  W )  ~  h‘ 


(6) 

(7) 


where  Re  =  Uc/v.  Differentiating  the  first  expression  with  respect  to  y  and  the 
second  with  respect  to  x  and  subtracting  eliminates  the  pressure  term.  Using  the 
definition  of  vorticity,  u>  =  vx  —  u„,  we  have, 


dw  du>  du>  _  \  ( d2u  d2u;\  df\  df2 

dt  +  U  dx  V  dy  Re  \  dx 2  dy 2  )  +  dy  dx 


(8) 


This  represents  the  modified  form  of  the  vorticity-transport  equation;  the  definition 
of  streamfunction  and  the  continuity  equation  remain  unaltered.  The  expressions 
for  /i  and  /2  are  obtained  by  examining  the  component  form  of  Eq.  (1): 


fi  =  R\  —  2vQ  —  yfi  —  xtt2, 

(9) 

/2  =  f?2  +  2uD  +  xil  —  yfi2, 

(10) 

where  R  =  Riii  +  R 2e2,  and  f  =  xej  -f  ye2. 

The  streamfunction,  'P,  is  defined  by  the  relations 

Vy  =  u  =  u  *  e x, 

(11) 

CN 

<*) 

II 

s> 

II 

n 

* 

1 

(12) 

A  disturbance  streamfunction  is  also  defined  [9]  such  that 

=  V(x,y)-y , 

(13) 

y  —  y  cos  q  —  x  sin  a, 

(14) 
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Figure  3.  Airfoil  with  Linear  and  Torsional  Springs  Affixed 

establishing 


=  COS  G  -f  1pv, 

(15) 

=  sin  a  —  ipx. 

(16) 

Ancillary  to  the  vorticity-transport  equation  is  the  Poisson  equation  for  ip  [8] 

V20  =  =  -<*>•  (17) 

The  airfoil  with  the  linear  and  torsional  springs  affixed  is  shown  in  Figure  3. 
The  governing  equations  for  the  airfoil  with  two  degrees  of  freedom,  expressed  in 
dimensional  form,  are 

m0h  +  Saa  +  Dhh  +  Kkh  =  Qh, 

Sah  +  Ia&  +  Daa  +  Ka{ot  —  Go)  =  Qa. 
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(18) 

(19) 


where  h  is  the  vertical  displacement,  m0  is  the  airfoil  mass  per  unit  span,  Sa  is  the 
mass  static  moment,  Ia  is  the  mass  moment  of  inertia,  Kh  and  Ka  are  the  spring 
stiffness  coefficients,  and  Dk  and  Da  are  the  structural  damping  coefficients.  Qk 
is  the  net  applied  aerodynamic  force  in  the  vertical  direction,  while  Qa  is  the  net 
applied  aerodynamic  pitching  moment. 

The  transformation  of  the  airfoil  equations  of  motion  to  a  first-order,  coupled 
system  proceeds  via  the  straightforward  substitutions,  yx  =  h,  y2  =  A,  y3  =  a,  and 
3/4  =  a,  leading  to 

Vi  =  J/2,  (20) 

2/2  =  ( Qh  -  say4  -  Dhy2  -  Khyi)/m0 ,  (21) 

2/3  =  !M,  (22) 

2/4  =  (Qa  ~  Say2  -  Da y4  -  Ka(y3  -  a0))/Ia.  (23) 

The  equations,  in  this  form,  are  then  integrated  in  time  using  a  fourth-order 
Runge-Kutta  algorithm  [33].  The  integration  takes  place  immediately  after  the  cal¬ 
culation  of  the  aerodynamic  coefficients  in  the  main  algorithm  and  proceeds  from 
time  t  to  t  +  At.  The  linear  and  angular  accelerations  obtained  are  then  returned  to 
provide  updated  values  for  the  aerodynamic  calculations. 

2.2  Coordinate  Transformation 

The  equations  of  motion  as  presented  are  applicable  to  a  Cartesian  coordinate 
system.  To  transform  the  equations  so  that  they  apply  to  a  general  curvilinear 
coordinate  system,  the  chain  rule  of  differentiation  is  employed.  The  transformation 
is  from  the  physical  domain  ( x,y )  to  the  computational  domain  (£,  17),  where 

Z  =  t(x,y),  y  =  v(x,y)-  (24) 
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The  direct  application  of  the  chain  rule  then  produces  the  derivative  operators 

(25) 

(26) 

Here,  the  subscripts  denote  differentiation  with  respect  to  time  or  spatial  coordinate. 
The  Jacobian  of  the  inverse  transformation  is  assumed  to  be  non- vanishing  [34]: 

J  =  x4y„  -  xvy(  ±  0.  (27) 

The  coordinate  transformation  results  in  the  following  form  for  the  Laplace  operator 

[8]: 

j2VV  =  Mil  +  Mrm  -  2^3 +  Ml  +  Mr,  =  (28) 

where  the  transformation  coefficients  are  functions  of  the  metrics  and  their  deriva¬ 
tives, 

4>\  =  Xi  +  y}, 

<t>7  =  x2v  +  y*, 

<h  =  Xii^  +  yty,,, 

<t> 4  =  J~l[MVnnxn  ~  *mVn)  +  Mmxn  ~  *&Vn)  +  2MxtvVn  ~  Vtn*v)]> 

=  -j~l{Myvr,Xi-xr)Vyt)  +  <j>2(yiiXi-xxyi)  +  2<h(xir,yi-yiTlxi)}. 

The  transformed  Laplacian  is  applied  to  the  Poisson  equation  for  the  stream- 
function  and  to  the  vorticity  transport  equation,  where  the  latter  equation  takes  the 
following  form: 

l  #>  A  A 

u ;,  +  uu>4  +  vuv  =  +  /i  -  /2, 


d__  d_  d_ 

dx~^d(+T,tdv' 
d_  _  d_  d_ 

dy  ~  ^vd(  +r?v<V 


(29) 


with 


“  =  y  =  j(y*  +  V>i),  *>  =  -y-  =  -j(y«  +  V><), 

/i  =  &/i*  +  Vvfivi  h  =  ^xfii  +  Vzhri- 

2.3  Aerodynamic  Coefficients 

The  streamfunction-vorticity  formulation  of  the  governing  equations  eliminates 
pressure  as  a  flow  variable.  Therefore,  an  alternate  means  of  calculating  pressure  is 
required  to  obtain  the  aerodynamic  coefficients.  The  procedure  employed  is  similar 
to  that  utilized  by  Mehta  [35].  The  pressure  is  integrated  along  a  path  normal  to 
the  airfoil  surface  to  establish  the  pressure  at  the  leading  edge;  the  integration  then 
proceeds  along  the  airfoil  upper  and  lower  surfaces  to  establish  the  pressure  at  these 
locations.  It  should  be  noted  that  the  pressure  at  the  leading  edge  is  not  required  for 
the  calculation  of  lift,  drag,  or  moment.  However,  in  other  applications  the  surface 
pressure  distribution  was  required,  so  the  present  technique  is  maintained.  The  basic 
approach  is  to  manipulate  the  components  of  the  momentum  equation  to  yield  p{  and 
pv.  The  components  of  the  momentum  equation  require  modification  to  account  for 
the  non-inertial  character  of  the  airfoil-fixed  coordinate  frame.  The  transformation 
proceeds  from  Eqs.  (6)  and  (7),  and  results  in 

«t  +  j(*„w€  -  +  j(wy«P»»)  =  ^v2«-/i,  (30) 

Vt  +  -  *<«„)  +  i (x*p„  -  xvPi)  =  -  /2,  (31) 

where  the  Laplacian  operator  now  refers  to  the  transformed  domain  as  in  Eq.  (28). 

Using  the  identities 


**VJu  +  y{V2u  =  +  (fow*), 

(32) 

x„V2u  +  y„V2u  =  j(-<t>3u>v  +  <f>  2U>f), 

(33) 
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Eqs.  (30)  and  (31)  are  manipulated  to  produce 

Pi  =  J  ~  +  *<(*«“*  +  y(Vv)  -  #„(x«u£  +  y£t>£)J 

~xih  -  y«/2  -  -  ytvu  (34) 

Pv  =  j  _  ^aW»)  +  +  yr>vv )  -  ’Mx»>“*  +  yn”c)] 

xvfi  yr,h  —  y^vt-  (35) 

On  the  airfoil  surface,  Eq.  (34)  simplifies  to 

Pi  ~  J^(&»w£  “  four,)  -  X{fi  -  y£/2.  (36) 

The  pressure  is  first  calculated  at  a  node  on  the  far-field  boundary  using  Bernoulli's 
equation, 

Poo  =  ^(1  ~  «L  -  *£>)■  (37) 

This  is  followed  by  itegration  along  a  line  of  constant  £  to  obtain  the  pressure  at  the 
leading  edge  of  the  airfoil, 


The  pressure  along  the  upper  and  lower  surfaces  is  then  computed  by  integration  of 
Eq.  (36), 

P(6)  =  /  hp(dt+p(Za).  (39) 

The  aerodynamic  coefficients  may  then  be  determined.  The  chord  and  normal 
forces  are  calculated  first,  from  which  the  lift  and  drag  coefficients  are  obtained 


Cc 

Cn 


2  £  pdy  -  ~  j>  udx, 
-2  fpdx-  ^-£udy. 


(40) 

(41) 
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The  lift  and  drag  coefficients  are  obtained  by  transforming  Cc  and  Cn  to  the  wind 
reference  frame 

Ci  =  C„  cos  a  —  Cc  sin  a, 

Cd  =  Cc  cos  a  +  Cn  sin  a. 

Finally,  the  pitching  moment  about  the  center  of  gravity  is  determined  via 

Cm  =  j>  rA  x  fthdA.  (44) 

where  r*x  is  the  vector  from  the  center  of  gravity  to  the  differentia]  area,  dA,  n  is 
a  unit  outward  normal  to  dA  ,  and  ft  is  the  total  force  (per  unit  area)  acting  on 
element  dA.. 

2-4  Boundary  Conditions 

The  boundary  conditions  on  the  surface  of  the  airfoil  are  given  by  Beran  and 
Lutton  [9],  where  is  specified  to  vanish 

$  =  ip  +  y  —  0,  (airfoil  surface).  (45) 

The  no-penetration  condition  (45)  and  the  no-slip  condition  require  that  and  \P„ 
vanish  on  the  surface,  which  is  equivalent  to  the  specification 

=  0,  (airfoil  surface).  (46) 

This  specification  is  used  to  develop  an  implicit  condition  for  surface  vorticity  from 
the  evaluation  of  the  streamfunction  at  the  surface  [9] 

4xVm+uj3  =  0.  (47) 


(42) 

(43) 
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The  specification  of  the  conditions  on  the  outer  computational  boundary  is 
complicated  by  the  fact  that  they  are  dependent  upon  12(f).  The  absolute  velocity 
in  the  far-field  is  =  cos  ae.\  +  sin  ae2  expressed  in  the  a-frame  coordinate  system. 
From  the  previous  analysis 


Voo  =  ~rr  +  “oo  +  1)  x  r, 
at 


so  the  relative  velocity  at  the  outer  boundary  is  then 

*»  = 

-  u&ex  +  Voo^.  (49) 

Examining  the  components  of  the  previous  equation  establishes  the  expression  for 
the  velocity  components  on  the  outer  computational  boundary 

Hoc  =  (cos  a  —  Ri  +  yO)ei  4-  (sin  a  —  R2  —  xl))e2.  (50) 

The  expressions  for  the  vorticity  and  streamfunction  on  the  outer  boundary 
are  then  obtained  by  differentiation  and  integration,  respectively, 

Woo  =  -21),  (51) 

*oo  =  (cos  a-  Ri)y  -  (sin  a  -  R2)x  +  ^l)(x2  +  y2).  (52) 

The  disturbance  streamfunction  on  the  outer  computational  domain  is  then 


V’oo  =  $00  -  y  =  -R\y  +  R2X  +  -l)(x2  +  y2). 
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2.5  Numerical  Implementation 

The  computational  coordinates,  £  and  r?,  axe  discretized  and  referenced  by  the 
indices  i  and  j,  respectively: 

1  <  i  <  /,  1  <  j  <  J. 

Spatial  derivatives  in  the  Laplacian  are  discretized  using  second-order  accurate, 
central-difference  expressions.  The  convective  terms  in  the  vorticity-transport  equa¬ 
tion  are  treated  with  a  hybrid  approach  [8].  The  convective  term  is  approximated 
using  central  differences,  while  a  second-order  accurate,  upwind  difference  is  applied 
to  the  other  convective  term,  uu>{.  This  approach  improves  the  smoothness  of  the 
vorticity  field  in  the  absence  of  explicit  artificial  viscosity. 

The  aerodynamic  equations  of  motion,  (Eqs.  (28)  and  (29)),  are  integrated 
in  time  using  a  semi-implicit  procedure  applying  LU  decomposition,  as  developed 
by  Beran  [8].  The  Poisson  equation  and  the  viscous  terms  in  the  vorticity-transport 
equation  are  treated  implicitly.  The  convective  terms  in  the  vorticity-transport  equa¬ 
tion  are  treated  explicitly  via  a  first-order  accurate  Euler  scheme. 

The  time-integration  procedure  is  written  in  delta  form,  with  the  correction 
vectors,  A "0  and  Anu;  defined  by 

0n+1  =  0n  +  (0n+1  -  0")  =  0n  +  An0, 

wn+1  =  wn  +  (wn+1  -  un)  =  wn  +  Anu>, 

where  n  is  the  temporal  index  such  that  t  =  nAt. 

The  Poisson  equation  relating  streamfunction  and  vorticity  (Eq.  28),  at  points 
away  from  the  airfoil  surface,  may  be  recast  using  operator  notation  as 

L0n+1  +  wn+1  =  0, 
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LAnip  +  A”u>  =  -(Ltpn  +  un)  =  Ni. 

(54) 

Decomposing  the  Laplace  operator  into  two  linear  operators,  L\  acting  on  the  interior 

nodes,  and  L2  acting  on  the  nodes  adjacent  to  the  airfoil  surface,  yields 

LxAnipi  +  L2AnrJ>t  +  Anu;j  =  Ni, 

(55) 

where  the  subscripts  i  and  s  refer  to  interior  and  surface  nodes,  respectively.  At  the 

nodes  on  the  airfoil  surface,  $  vanishes,  leading  to  the  discrete  form  of  the  boundary 

condition 

0"+1  +  y  =  0, 

A "0.  =  -«?  +  »)  3  AT,. 

(56) 

Likewise,  the  discrete  boundary  condition  for  vorticity  on  the  airfoil  surface 

expressed  as 

may  be 

L30,n+1+u'r1=°i 

£3A"0i  +  Anu;,  =  —  (L30”  +  v?)  =  N3. 

(57) 

The  vorticity-transport  equation,  evaluated  at  nodes  away  from  the  surface,  is  writ- 

ten  in  discrete  form  as 

Anu,  -  ~Lun+l  =  -(uw{  +  vvv)nA t  +  (A  -  f2)At, 
tit 

or  when  reexpressed  in  delta  form 

At  At 

Anu>  —  -—LAnu>  =  —(itujt  4-  vu>v)nAt  +  — Lu/1  -f  (/1  -  f2)At  =  iV4. 
ne  nt 

The  operator  G  is  defined  such  that 

GA”W  .  [/  -  f  l]  AV 
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GiAnUi  +  Gt  A"u>.  =  N4.  (58) 

Provided  that  At  is  held  constant,  Gi  and  Gt  «j*e  constant  matrices.  The  resulting 
blocked  system  of  equations,  (55-57)  .  can  be  reduced  to  a  single,  banded  system  for 
Anr/>,- 

-4An0,  =  N4  +  GiihN*  -  Nt)  -  G.N3 ,  (59) 

where 

A  =  -(G.Ii  +  G,L3). 

The  matrix  A  is  time-invariant  (for  a  constant  time  step  and  Reynolds  number),  and 
is  therefore  decomposed  into  the  product  of  lower  and  upper  triangular  matrices  at 
the  start  of  the  time-integration  procedure.  Once  and  0n+1  are  calculated,  the 
vorticity  can  be  updated  by  applying  Eqs.  (55)  and  (57): 


Anu>,  =  Ni-  IiAnV>,  -  L2N2,  (60) 

Anu>.  =  N3  -  IsA’Vi.  (61) 

The  far-field  boundary  conditions,  implemented  in  delta  form,  are 

Anu>oo  =  -2Anfi,  (62) 

An0oo  +  AwMoo  =  -AnRiy  +  AnR2x  +  ^Anfi(x2  +  y2  -  4).  (63) 


Finally,  the  velocities  and  aerodynamic  coefficients  are  calculated,  and  the  integra¬ 
tion  of  the  airfoil  equations  of  motion  accomplished.  The  aerodynamic  coefficients 
are  calculated  immediately  after  the  resolution  of  xl>  and  u>  and  are  treated  as  con¬ 
stants  over  the  interval  t  to  t  +  At.  The  integration  of  the  airfoil  equations  of 
motion  (18,19)  proceeds  over  the  same  interval  by  subdividing  this  interval  to  apply 
the  Runge-Kutta  integration,  using  the  constant  values  for  the  applied  aerodynamic 
loads.  The  resulting  linear  and  angular  accelerations  are  then  returned  to  provide  up- 
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Figure  4.  Flowchart  Detailing  Implementation  of  Aerodynamic  and  Structural 
Models 

dated  values  for  the  aerodynamic  calculations  at  the  beginning  of  the  next  timestep. 
This  introduces  a  phase  lag  of  up  to  At  in  the  application  of  the  loads,  however  this 
is  consistent  with  the  first-order  accuracy  in  time  of  the  overall  scheme.  The  imple¬ 
mentation  of  the  structural  model  is  detailed  in  Figure  4.  This  procedure  has  been 
established  as  robust  by  comparing  airfoil  response  to  a  prescribed  forcing  function. 
The  validation  procedure  is  documented  in  Appendix  A. 
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III.  Results  and  Conclusions:  Part  I 


S.l  Grid  Refinement  and  Sensitivity 

Numerical  experiments  were  conducted  with  a  NACA  0012  airfoil  using  a  vari¬ 
ety  of  C-grids.  All  grids  were  constructed  with  a  hyperbolic  grid  generation  algorithm 
developed  by  Barth  and  Kinsey  [36].  The  grid  generation  algorithm  proceeds  from 
a  point  distribution  specified  on  the  airfoil  surface.  Hyperbolic  partial  differential 
equations  are  employed  to  establish  the  normal  point  distribution,  providing  a  nearly 
orthogonal  mesh.  Grids  1  and  2  are  shown  in  Figures  5  and  6.  Table  1  provides  the 
numerical  characteristics  of  the  grids,  displaying  the  grid  dimensions,  the  number  of 
nodes  on  the  upper  surface  (Nu),  lower  surface  (jV;),  and  wake  (Nw),  the  tangential 
node  spacing  at  the  leading  edge  (As/e),  and  trailing  edge  (Aste),  and  node  spacing 
normal  to  the  surface  (As„).  Also  provided  in  Table  1  is  the  distance  to  the  outer 
boundary  ( Rj ),  in  chord  lengths.  Because  of  computational  constraints,  the  majority 
of  numerical  investigations  were  performed  using  Grid  1;  selected  data  points  were 
reexamined  using  other  grids. 

To  assess  the  impact  of  grid  refinement,  the  Strouhal  number  and  the  critical 
Reynolds  number.  Rec-u,  defined  as  the  Reynolds  number  at  which  the  flow  becomes 
unsteady,  were  determined  for  the  fixed  airfoil  case.  The  Strouhal  number  was 
determined  by  measuring  the  period  of  the  oscillations  in  Ct  at  Re  =  1200  for  each 
grid.  Rtcrit  was  determined  by  starting  from  an  equilibrium  solution,  and  integrating 
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Table  1.  Grid  Characteristics 
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C-Grid 

IfXfl 

St 

R^crit 

1 

59x25 

0.93 

550 

ID 

99x40 

0.95 

760 

IKD 

139x50 

0.98 

810 

179x60 

0.99 

850 

Table  2.  Strouhal  no.  at  Re  =  1200  and  Hopf  points,  a  =  0.2  rad 

in  time  to  assess  the  stability  of  the  initial  solution,  as  in  Beran  and  Lutton  [9],  The 
results  are  displayed  in  Table  2.  The  Strouhal  number  shows  a  clear  convergence, 
varying  by  slightly  more  than  one  percent  between  grids  3  and  4.  The  /?eCTlt  data 
show  a  similar  convergence,  varying  by  less  than  five  percent  between  grids  3  and 
4.  This  increased  variation  is  expected,  as  the  determination  of  Re^u  is  less  precise 
than  the  Strouhal  number  determination. 

In  addition,  the  effect  of  varying  the  placement  of  the  outer  boundary,  Rj, 
was  explored.  This  was  accomplished  by  constructing  a  series  of  grids,  each  with 
parameters  identical  to  those  of  Grid  2,  except  for  the  variation  of  Rd.  The  Strouhal 
number  is  determined  in  each  case,  again  at  Re  =  1200.  The  results  are  displayed 
in  Figure  7.  The  Strouhal  number  shows  very  little  variation  for  Rd  >  4.  It  should 
be  noted  that  grid  parameters  play  an  important  role  in  this  type  of  investigation. 
For  instance,  if  the  trailing  edge  spacing,  Aste,  is  increased  to  .005  for  Grid  2,  the 
Strouhal  number  decreases  to  0.41,  clearly  inconsistent  with  the  stated  results.  This 
is  attributed  to  the  failure  to  resolve  flow  structure,  namely  vortices,  in  the  vicinity 
of  the  trailing  edge.  At  increased  domain  sizes  (Rd  >  10)  a  slight  decrease  in  the 
Strouhal  number  is  noted,  apparently  due  to  an  insufficient  number  of  wake  points. 
Considering  these  effects,  it  is  a  necessity  to  perform  at  least  a  limited  grid  refinement 
study  when  addressing  similar  topics. 
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3.2  Variation  of  Natural  Frequencies:  Single  and  Dual  Axes 

The  procedure  employed  to  determine  the  approximate  location  of  Hopf  points 
is  explained  here.  Steady  solutions  were  obtained,  using  time  integration,  at  very 
low  Reynolds  numbers  (Re  <  100).  Continuation  in  Reynolds  number  was  then 
employed,  i.e.,  the  Reynolds  number  was  gradually  increased  in  increments  between 
10  and  100,  depending  on  the  rate  of  convergence  of  the  previous  solution,  until  an 
oscillatory  solution  was  obtained.  The  instantaneous  change  in  Reynolds  number 
introduces  perturbations  into  the  numerical  flowfield,  which  are  then  either  numer¬ 
ically  damped,  or  evolve  to  an  unsteady  flow.  The  value  of  Reynolds  number  at 
which  the  flow  becomes  unsteady  is  denoted  as  Rea-n,  representing  the  approximate 
location  of  the  Hopf  point.  This  procedure  is  not  exact;  the  main  interest  being  the 
shift  in  the  Reynolds  number  between  the  fixed  and  moving  airfoils.  The  principal 
difficulty  in  this  type  of  investigation  is  the  amount  of  computational  time  required 
to  accurately  determine  the  Hopf  point.  As  the  structural  and  inertial  properties  of 
the  airfoil  are  varied,  the  described  procedure  must  be  reapplied.  However,  as  the 
Hopf  point  is  approached,  the  aerodynamic  damping  becomes  less  effective,  and  a 
large  number  of  time  steps  axe  required  to  assess  the  stability  of  the  computation. 

Numerical  investigations  were  performed  for  a0  =  0.1  rad  and  oD  =  0.2  rad, 
where  a0  is  the  value  of  a  at  which  the  torsional  spring  delivers  no  moment  (cf.  Eq. 
(19)).  For  each  of  these,  a  variety  of  cases  were  considered.  The  natural  frequency 
in  pitch,  uja  =  yjKa/Ia ,  was  varied  from  5  to  50;  similarly,  the  variation  of  the 
natural  frequency  in  plunge,  =  \jKhlm0,  was  varied  across  the  same  range. 
These  investigations  are  undertaken  with  both  axes  active  (dual  axes),  and  with  the 
complimentary  axis  disabled  (single  axis),  to  assess  the  influence  of  the  the  individual 
axes,  and  the  axes  in  tandem,  on  the  location  of  the  Hopf  point.  In  the  dual  axes 
case,  u>a  was  varied,  while  was  held  constant.  The  structural  coupling  term,  S0, 
was  initially  set  to  zero. 
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Figure  5-  Grid  1  (59x25) 


The  results  for  the  a0  =  0.2  case  are  shown  in  Figure  8.  The  onset  of  peri¬ 
odic  motion  for  the  fixed  airfoil,  determined  by  the  procedure  outlined  above,  was 
determined  to  occur  at  approximately  Re  =  475.  This  is  reasonably  close  to  the 
value  predicted  by  a  direct  method  for  calculation  of  the  Hopf  point,  which  yields 
Re  =  550  for  the  fixed  airfoil  [9j.  The  direct  method  involves  solving  an  extended 
set  of  equations,  augmented  to  capture  the  eigenvalue  migration  as  the  Reynolds 
number  changes  (cf.  Appendix  C).  In  the  direct  method,  the  Reynolds  number 
represents  an  unknown,  and  the  solution  delivers  the  Reynolds  number  for  which  a 
pair  of  complex  conjugate  eigenvalues  are  migrating  across  the  imaginary  axis  (cf. 
Appendix  C).  This  represents  a  mathematical  description  of  Hopf  bifurcation.  Two 
methods  are  therefore  available  to  determine  the  Hopf  bifurcation  point:  the  direct 
method  and  time  integration.  See  Beran  and  Lutton  [9]  and  Jackson  [3]  for  a  detailed 
presentation  of  the  direct  method.  The  difference  in  the  two  methods  suggests  that 
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Figure  6.  Grid  2  (99x40) 

the  structure  of  the  Hopf  bifurcation  is  subcritical ,  a  question  explored  further  in 
Sections  3.3  and  3.4.  With  both  axes  active,  there  was  a  significant  reduction  in  the 
Reynolds  number  at  which  the  flow  became  periodic,  to  approximately  Re  =  300, 
across  a  broad  range  of  uia.  The  exception  occurs  for  values  of  uQ  below  10.  This  ef¬ 
fect  is  attributed  to  the  smaller  value  of  Ka  employed.  The  lower  value  of  Ka  allows 
the  airfoil  to  equilibrate  to  a  new,  lower  angle  of  attack,  for  which  a  larger  value  of 
the  critical  Reynolds  number  is  expected.  This  is  consistent  with  earlier  results  [9]. 
The  single-axis  experiments  are  also  displayed  in  Figure  8.  The  most  notable  result 
is  that,  with  the  pitch  axis  fixed,  the  effect  of  the  vertical  axis  on  the  Hopf  point 
location  is  negligible  across  the  range  of  ufo,  except  at  uv,  =  5,  where  a  reduction 
to  Re  =  350  is  noted.  This  result  is  due  to  resonance,  since  the  angular  frequency 
of  the  aerodynamics,  determined  by  examining  the  fixed  airfoil  at  Re  =  500  and 
a  =  0.2  rad,  is  5.06  rad/sec.  It  should  be  noted  that  the  angle  of  attack  was  fixed 
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Figure  7.  Variation  of  Strouhal  no.  with  Domain  Size 

in  this  case  (at  a  —  0.2),  unlike  the  previous  case,  where  a  was  not  fixed.  The 
results  of  a  second,  single-axis  experiment  are  also  shown  in  Figure  8.  Here,  the 
vertical  axis  was  fixed,  and  the  pitch  axis  was  active.  A  reduction  in  the  critical 
Reynolds  number  was  again  noted,  with  the  same  trend  observed  for  u>h  less  than 
10.  The  differences  in  the  critical  Reynolds  number  noted  in  the  dual  and  single-axis 
experiments  indicate  that  the  axes  become  coupled  via  the  aerodynamics  even  in  the 
absence  of  direct  structural  coupling.  The  aerodynamic  coupling  occurs,  and  has  an 
effect  on  the  onset  of  periodic  motion,  despite  the  fact  that  the  vertical  axis,  alone, 
has  a  minimal  impact  for  >  10. 

The  same  set  of  numerical  experiments  was  then  repeated  at  a0  =  0.1  rad  to 
determine  if  the  same  trends  are  evident  at  a  lower  angle  of  attack.  See  Figure  9.  The 
critical  Reynolds  number  for  the  fixed  airfoil  occurs  at  Re  =  1600  when  applying  Re 
continuation,  while  the  direct  solver  yields  Re  =  2320.  The  increased  difference  be¬ 
tween  these  two  numbers,  when  compared  with  the  a  =  0.2  case,  again  indicates  the 
subcritical  nature  of  the  bifurcation,  and  that  the  bistable  region  grows  as  the  angle 
of  attack  decreases.  This  may  explain  the  results  observed  for  the  vertical  axis  case: 
no  effect  upon  the  location  of  the  critical  Reynolds  number  for  =  50,  however,  at 
lower  values  of  u the  value  of  the  critical  Reynolds  number  is  actually  increased. 
The  differences  in  the  critical  Reynolds  number  observed  between  the  fixed  airfoil 
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Figure  8.  Critical  Reynolds  no.  vs  w,  a  =  0.2  rad 

and  vertical-axis  case  indicate  that  the  effect  of  the  structural  model  is  the  migra¬ 
tion  of  the  stable  and  unstable  periodic  orbits  associated  with  the  bistable  region. 
That  is,  the  location  of  the  periodic  orbits  (both  stable  and  unstable  branches)  in 
the  bistable  region  are  more  sensitive  to  the  addition  of  a  structural  model  than 
is  the  actual  location  of  the  Hopf  point.  This  may  occur  because  these  orbits  are 
inherently  periodic,  and  hence  immediately  have  the  potential  to  excite  the  airfoil. 
See  Figure  15  for  an  indication  of  this  effect.  Since  analysis  reveals  the  vertical  axis, 
singularly,  should  have  no  direct  influence  upon  the  location  of  the  Hopf  bifurcation 
point,  the  numerical  results  presented  are  attributed  to  the  migration  of  periodic 
orbits  in  the  bistable  region. 

The  dual-axis  experiments  show  the  same  trend  (reduction  in  the  critical 
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Figure  9.  Critical  Reynolds  no.  vs  u,  a  =  0.1  rad 

Reynolds  number)  as  the  a0  =  0.2  case,  except  that  the  reduction  is  of  greater 
magnitude.  It  is  observed  that  the  initial  onset  of  oscillations  begins  at  the  same 
Reynolds  number,  Re  =  300,  for  both  the  a  =  0.1  and  a  =  0.2  cases  with  u>a  >  10. 
The  same  trend  at  lower  values  of  ua  is  evident  as  the  airfoil  equilibrates  to  a  lower 
angle  of  attack.  The  differences  in  the  value  of  critical  Reynolds  number  between 
the  dual-axis  results  and  the  pitch-axis  results  are  again  apparent,  but  are  also  of 
greater  magnitude  than  in  the  a0  =  0.2  rad  case. 

The  effect  of  structural  coupling,  Sa,  was  also  examined.  The  value  for  the 
previously  discussed  data  was  Sa  =  0.0.  The  value  of  Sa  was  increased  to  0.1,  and 
again  to  0.2.  With  both  axes  active,  =  50  and  =  10,  there  was  no  significant 
change  in  the  critical  Reynolds  number  for  either  value  of  Sa.  Since  the  single-axis 
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and  dual-axis  experiments  demonstrated  the  influence  of  coupling,  it  is  probable  that 
the  axes  are  strongly  coupled  aerodynamically,  regardless  of  the  presence  of  direct 
structural  coupling. 

3.3  Influence  of  Structural  Damping 

The  influence  of  structural  damping,  Da  and  Dh,  was  explored  next.  Since  the 
onset  of  oscillatory  flow  is  closely  related  to  the  ratio  of  the  convective  and  dissipative 
phenomena  inherent  in  the  flow  via  the  Reynolds  number,  it  is  of  interest  to  deter¬ 
mine  the  effect  on  the  onset  of  oscillations  if  an  additional  dissipative  mechanism  is 
applied.  The  results  are  shown  in  Figures  10,  11  and  12.  These  figures  present  the 
aerodynamic  response,  as  indicated  by  Q,  and  the  structural  response  of  the  airfoil, 
indicated  by  a  and  h,  as  the  levels  of  structural  dissipation  are  varied.  A  Fourier 
analysis  [37]  is  also  presented  to  quantify  the  frequency  content  of  the  aerodynamic 
response.  The  frequency  is  presented  as  Hertz  (Hz),  here  measured  as  cycles  per 
aerodynamic  time  unit.  The  amplitude  is  scaled  such  that  a  constant  signal  will 
provide  an  amplitude  equal  to  its  value,  while  a  single-frequency,  harmonic  signal 
will  provide  an  amplitude  equal  to  half  of  its  peak-to-peak  amplitude. 

The  results  indicate  that  the  effect  of  increased  levels  of  structural  dissipation 
are  dependent  upon  the  Reynolds  number.  The  damping  factor,  (,  is  defined  such 
that  £  =  D/Dcrit,  where  Do-a  represents  the  critical  damping  level,  here  equal  to 
2\/m0Kh  for  the  vertical  axis  or  2y/ IaKa  for  the  pitch  axis.  All  numerical  experi¬ 
ments  were  performed  with  the  damping  factor  equal  for  both  axes,  hence  a  single 
value  of  (  is  given  for  each  result  rather  than  stating  a  value  for  each  axis.  Fig¬ 
ure  10  displays  the  results  at  F.e  =  450,  (  =  0.5,  representing  an  underdamped  case. 
The  unsteady  flow  results  in  oscillatory  airfoil  motion.  As  the  damping  levels  are 
increased  to  C  =  2.0  for  the  same  Reynolds  number  (Figure  11),  the  airfoil  motion 
is  damped  and  eventually  a  steady  state  is  achieved.  The  behavior  is  different  at 
He.  =  650,  however.  When  the  damping  levels  were  increased  to  (  =  2.0,  the  flow 
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field  and  the  airfoil  remain  oscillatory,  as  shown  in  Figure  12.  The  difference  in  these 
two  cases  appears  to  be  related  to  the  Hopf  bifurcation  structure.  In  the  first  case, 
Re  =  450,  the  Reynolds  number  is  located  in  the  bistable  region  (See  Figure  15), 
and  the  added  dissipation  has  the  effect  of  changing  branches  from  the  (stable)  limit 
cycle  solution  to  the  (stable)  equilibrium  solution.  Conversely,  Re  =  650  is  beyond 
the  bistable  region  and  the  added  dissipation  cannot  damp  the  airfoil  motion,  as  the 
equilibrium  solution  in  this  case  is  unstable.  Further  numerical  experiments  were 
performed  with  £  ranging  up  to  50.  Even  at  these  high  levels,  the  results  remained 
the  same:  at  Re  =  450  an  equilibrium  solution  was  obtained,  while  at  Re  =  650 
oscillations  persisted  in  C/,  a,  and  h. 

Once  the  Reynolds  number  has  increased  past  the  upper  bound  of  the  bistable 
region,  the  additional  dissipation  is  seen  to  reduce  the  amplitude  of  the  airfoil  os¬ 
cillations,  but  does  not  eliminate  them,  even  at  very  high  damping  levels  (£  =  50). 
This  result  is  explained  by  the  nature  of  the  additional  dissipation  applied.  That  is, 
viscosity  in  the  flowfield  has  an  impact  whenever  velocity  gradients  develop  in  the 
flow,  principally  in  the  boundary  layer  and  the  wake.  In  contrast,  the  structural  dis¬ 
sipation  is  only  effective  when  the  airfoil  is  moving,  and  is  proportional  to  a  and  h. 
Therefore,  the  impact  of  the  structural  dissipation  only  becomes  apparent  after  the 
airfoil  begins  to  oscillate  in  either  pitch  or  vertical  displacement.  In  this  Reynolds 
number  range  (650),  the  lift  coefficient,  for  example,  will  remain  oscillatory  for  any 
level  of  damping  (i.e.  the  fixed  airfoil  exhibits  oscillatory  lift). 

3.4  Impact  of  the  Structural  Model  on  the  Solution  Space 

In  this  section,  comparisons  are  made  between  the  solution  spaces  for  the  cases 
of  the  fixed  and  the  moving  airfoils.  The  general  term  “solution  space”  is  employed 
because  it  refers  to  a  variety  of  results:  vorticity  in  the  near  wake,  bifurcation  struc¬ 
ture,  and  Ci  phase  plots.  In  general,  the  addition  of  structural  elements  to  the  airfoil 
tends  to  produce  a  more  dynamically  varying  result:  additional  frequency  content, 
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migration  of  bifurcation  points,  and  an  inherently  different  behavior. 

The  impact  of  the  structural  model  on  the  flowfield  is  shown  in  Figures  13 
and  14.  Both  solutions  were  obtained  by  integrating  from  the  same  initial  state  at 
a  =  0.2  rad  using  Grid  2.  In  the  case  of  the  fixed  airfoil,  a  steady  flow  ensues. 
A  small,  stationary  vortex  is  seen  at  the  trailing  edge.  The  dissipative,  convective 
and  generative  forces  are  in  balance:  the  vorticity  generated  in  the  boundary  layer  is 
being  convected  to  the  stationary  vortex  at  the  same  rate  vorticity  is  being  dissipated 
from  the  same  structure.  In  contrast,  Figure  14  displays  the  vorticity  when  the  airfoil 
is  permitted  to  move.  The  aerodynamic  response  is  heightened,  and  vortex  shedding 
occurs.  The  airfoil  is  pitching  nose  down  at  an  instantaneous  rate  of  0.1  rad/sec, 
a  large  region  of  (negative)  vorticity  generated  in  the  boundary  layer  on  the  upper 
surface  has  convected  into  the  lower  wake.  Behind  this  vortex,  also  in  the  lower  wake, 
is  a  second  region  of  (positive)  vorticity  initially  generated  in  the  boundary  layer  on 
the  lower  surface.  The  differences  in  the  vorticity  patterns  clearly  demonstrate  that 
the  movement  of  the  airfoil  induces  the  development  of  unsteady  flow  before  that 
which  would  occur  for  the  fixed  airfoil. 

As  discussed  earlier,  the  structure  of  the  Hopf  point,  when  examined  using 
the  procedure  outlined,  appears  subcritical.  The  perturbations  introduced  into  the 
flowfield,  when  continuation  in  Reynolds  number  is  applied,  appear  to  initiate  the 
appearance  of  undamped,  limit-cycle  oscillations  at  a  lower  Reynolds  number  them 
would  otherwise  be  the  case.  In  contrast,  if  the  steady-state  solver  of  Beran  [8]  is 
applied,  and  an  equilibrium  solution  is  used  as  the  initial  state,  then  the  onset  of 
periodic  flow  is  delayed  with  respect  to  Reynolds  number.  Figure  15  illustrates  the 
structure  of  the  subcritical  bifurcation.  The  results  were  obtained  using  Grid  1  at 
a  =  0.2  rad.  The  solid  line  represents  stable  equilibrium  solutions,  obtained  using 
the  equilibrium  solver  of  Beran  [8].  The  equilibrium  solutions  become  unstable  at 
the  Hopf  point  (Re  =  550),  as  determined  by  the  direct  Hopf  solver,  and  the  solid 
line  is  replaced  by  the  dashed  line.  An  eigenvalue  analysis  at  Re  =  550  further 
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demonstrates  a  conjugate  pair  of  eigenvalues  migrating  across  the  imaginary  axis 
[9].  In  contrast,  no  Hopf  point  is  exhibited  along  the  path  of  equilibrium  solutions 
where  the  upper  and  lower  branches  of  the  periodic  orbits  appear  to  converge.  The 
outer  branches  of  Figure  15  represent  the  C/  extrema  of  the  limit-cycle  oscillations 
for  the  fixed  airfoil,  the  vertical  axis,  and  dual  axes  respectively.  The  solutions  along 
these  branches  are  obtained  by  continuation  from  a  previous,  oscillatory  solution. 
The  results  for  the  fi::ed  airfoil  and  the  vertical  axis  show  a  small  difference,  while 
the  branch  associated  with  the  dual-axes  case  shows  a  considerable  shift,  clearly 
indicating  the  impact  of  the  structural  model  on  the  location  of  the  stable,  periodic 
orbits.  It  should  be  noted  that,  given  the  subcritical  structure  of  the  bifurcation,  it 
becomes  problematic  to  determine  if  the  structural  model  has  an  effect  on  the  Hopf 
point,  or  whether  the  principal  effect  is  the  migration  of  the  path  of  periodic  orbits 
associated  with  the  bistable  region.  In  a  practical  setting,  the  question  may  not  be 
of  paramount  importance,  since  some  level  of  perturbation  is  likely  in  any  event: 
the  examination  of  the  periodic  orbits  in  the  bistable  region  would  then  be  of  more 
practical  importance. 

In  order  to  examine  the  dynamics  of  the  numerical  solutions  and  to  capture 
the  approach  to  attractors  in  the  solution  space  (i.e.,  equilibrium  or  fixed  points  and 
iimit  cycles),  phase  plots  are  introduced.  In  a  C;  phase  plot,  for  instance,  Ci(t)  is 
plotted  against  C/(<  +  St),  where  St  is  chosen  to  be  on  the  order  of  one-tenth  of  the 
period.  Thus  for  these  representations,  an  equilibrium  solution  is  represented  by  a 
single  point,  while  a  limit  cycle  represents  a  single,  closed  orbit.  The  shape  of  the 
orbit  is  dependent  upon  the  value  chosen  for  St.  Figure  16a  shows  the  evolution 
towards  a  limit  cycle  for  the  fixed  airfoil,  with  time  integration  proceeding  from  an 
equilibrium  solution  at  Re  =  1200,  a  =  0.2  rad,  using  Grid  2  (Rj  =  14).  The  phase 
plot  displays  a  trajectory  emanating  from  the  center  (equilibrium  state),  which  then 
overshoots  the  stable  orbit  before  contracting  back  to  a  fixed  amplitude.  Figure  16b 
indicates  the  corresponding  time  history.  Contrast  this  solution  with  the  case  of 
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the  airfoil  with  the  vertical  sods  active  (a;*,  =  6  rad/s),  displayed  in  Figure  17.  Time 
integration  proceeds  from  the  same  equilibrium  solution.  In  this  case,  the  progression 
to  a  set  of  bounding  orbits  is  less  regular,  and  the  solution  never  cleanly  approaches 
an  orbit  of  constant  amplitude.  This  is  also  indicated  in  the  time  history,  where 
additional  frequency  content  and  reduced  amplitude  is  evident  when  compared  to 
the  prior  case,  as  indicated  in  the  Fourier  analysis  of  Cj. 

It  is  ot  interest  to  consider  the  difference  in  the  behavior  of  the  two  cases  in 
the  neighborhood  of  the  equilibrium  solution,  i.e.,  to  examine  the  behavior  of  the 
small  perturbations  which  ultimately  bifurcate  the  solution  to  another  branch,  as  in 
Figure  15.  The  phase  plot  of  this  regime,  representing  trajectories  in  the  immediate 
vicinity  of  the  equilibrium  solution,  is  shown  in  Figure  18  for  the  fixed  airfoil.  The 
choppy  character  of  the  plot  is  due  to  the  finite  decimal  representation  of  C/,  but 
the  general  trend  remains  obvious.  The  trajectories  appear  to  construct  an  inner 
structure,  as  C;  is  slightly  perturbed  about  the  equilibrium  state.  However,  after 
leaving  this  inner  region,  the  trajectory  no  longer  intersects  itself  (until  much  later  as 
it  approaches  a  limiting  orbit),  and  proceeds  to  spiral  outward  in  a  regular  fashion. 
Contrast  this  with  the  result  for  the  case  with  the  vertical  axis  active,  shown  in 
Figure  19.  In  this  case,  the  trajectory  also  spirals  outward  from  the  neighborhood 
of  the  equilibrium  solution.  However,  the  trajectory  appears  to  double  back  and 
intersect  itself  repeatedly,  indicating  the  additional  frequency  content  of  the  solution. 
While  it  may  appear  macroscopically  that  the  difference  between  the  two  cases  are 
initially  negligible,  with  differences  appearing  only  over  extended  time  intervals, 
these  results  indicate  the  nature  of  the  bifurcations  are  fundamentally  different  and 
that  the  small  perturbations  about  an  equilibrium  solution  behave  differently. 

The  next  plot  (Figure  20)  was  extracted  from  the  same  data  set  used  to  con¬ 
struct  Figure  15  (Grid  1,  dual  axes,  u>a  =  25  rad/s,  =  6  rad/s).  It  is  essentially  a 
cross  section  of  Figure  15  and  is  intended  to  exhibit  the  evolution  towards  a  stable 
limit-cycle  solution  in  the  bistable  region  (Re  =  380).  Time  integration  proceeds 
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from  the  center  outward  in  the  phase  plot  (a);  the  corresponding  time  history  being 
shown  in  (b).  The  interesting  point  is  the  apparent  influence  of  the  unstable  limit- 
cycle  branch  which  connects  the  stable  limit-cycle  branch  to  the  Hopf  point.  A  limit 
cycle  of  lesser  amplitude  than  the  final,  stable  limit  cycle  appears  to  be  approached 
during  the  time  integration.  However,  since  such  a  solution  must  be  unstable,  the 
evolution  continues  to  the  larger  amplitude  limit  cycle.  It  is  noteworthy  that  this 
phenomenon  appears  uniquely  in  the  bistable  region,  and  was  not  observed  outside  of 
this  region,  see  for  example  Figure  16.  Behavior  of  this  type  has  also  been  observed 
in  flutter  experiments,  see  Yang  and  Zhao  [11]. 

The  final  plot  in  this  series,  shown  in  Figure  21,  is  provided  to  indicate  the 
range  of  dynamic  behavior.  The  plot  represents  a  case  with  both  axes  active,  at 
Re  =  950,  and  a  =  0.2,  with  uQ  =  =  25  rad/s,  utilizing  Grid  2.  The  phase 

plot  indicates  the  complex,  multiple  frequency  solutions  available  to  the  augmented 
dynamic  system.  The  multiple  frequencies  are  also  apparent  in  the  time  history  and 
in  the  corresponding  Fourier  analysis. 

3.5  Conclusions 

The  influence  of  the  structural  model  of  the  airfoil  on  the  location  of  the  Hopf 
bifurcation  structure  has  been  demonstrated.  In  general,  the  structural  model  allows 
the  development  of  sustained,  limit-cycle  oscillations  in  the  flowfield  to  develop  at 
a  lower  Reynolds  number  than  would  be  the  case  for  a  fixed  airfoil.  This  effect  is 
most  evident  when  the  structural  model  incorporates  multiple  degrees  of  freedom. 
Coupling  between  the  axes  can  occur  through  either  direct  structural  coupling,  or 
indirectly  via  the  aerodynamics.  Single-axis  experiments  demonstrate  that  the  ver¬ 
tical  axis  has  a  minimal  impact  on  the  value  of  Reynolds  number  at  which  the  flow 
becomes  oscillatory  in  the  absense  of  a  resonance  effect,  while  the  pitch  axis  has  a 
more  pronounced  effect.  The  effect  of  structural  damping  is  most  pronounced  in  the 
bistable  region,  beyond  this  region,  the  structural  damping  impacts  the  airfoil  motion 
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more  strongly  than  the  aerody manic  flowfield,  which  will  remain  unsteady  regardless 
of  the  structural  damping  level.  The  structure  of  the  bifurcation  is  subcritical,  as 
the  existence  of  perturbations  in  the  flowfield  can  initiate  sustained  oscillations  at 
a  Reynolds  number  lower  than  would  occur  when  starting  from  an  initial  equilib¬ 
rium  state.  The  influence  of  the  structural  model  may  be  more  pronounced  on  the 
branch  of  periodic  solutions  in  the  bistable  region  than  on  the  Hopf  point  itself.  In 
addition,  the  structural  model  is  seen  to  provoke  additional  frequency  content  in  the 
time  history  of  Cj,  for  instance.  The  addition  of  structural  elements  may  induce 
secondary  bifurcations  which  perturb  the  solution  path  from  a  single,  closed  orbit  in 
phase  space  to  a  more  complicated  orbit.  This  is  true  in  the  long-term  time  behavior 
and  in  the  immediate  evolution  away  from  an  equilibrium  solution. 

Grid  sensitivities  to  critical  parameters  such  as  the  Strouhal  number  and  crit¬ 
ical  Reynolds  number  are  considerable.  Relatively  minor  changes  in  grid  spacing 
or  domain  size  may  have  a  sizeable  effect  on  numerical  results.  At  a  minimum,  a 
limited  grid  refinement  study  is  therefore  required  to  establish  consistent  results. 
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Figure  15.  Hopf  Bifurcation  Structure:  a  =  0.2  rad,  Stable  Equilibrium  Path  (solid 
line),  Unstable  Equilibrium  Path  (dashed  line),  Dual  Axes  (O),  Fixed 
Axes  (©),  Vertical  Axis  (A) 
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Figure  19.  Vertical  Axis  Case,  Re  =  1600,  a  =  0.25  rad;  C/  Phase  Plot  Near 
Equilibrium 
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Figure  20.  Dual  Axes  Case,  Re  =  380,  a0  =  0.2  rad;  (a)  Ci  Phase  Plot  (b)  Time 
History  (c)  Fourier  Analysis47 


IV.  Analysis:  Part  II 


4-1  Turbulence  Model  Implementation 

One  of  the  most  difficult  and  challenging  problems  in  modern  computational 
fluid  dynamics  is  the  correct  modeling  of  turbulent  flow.  According  to  Bushnell 
[38]  the  importance  of  turbulence  modeling  is  so  paramount  that  it  represents  a 
primary  pacing  technology  in  the  development  of  CFD.  This  is  particularly  true  for 
complex  geometries  or  flowfields:  separated  and  recirculating  flows,  dynamic  flows, 
and  surfaces  with  curvature  or  comers  [14]. 

The  turbulence  model  employed  here  is  a  modified  version  of  the  algebraic, 
eddy-viscosity  model  of  Baldwin  and  Lomax  [13]  as  modified  by  Knight  and  Visbal 
[14].  The  primary  reason  for  selecting  the  Baldwin- Lomax  model  is  the  relatively 
low  computational  cost  and  ease  of  application.  In  this  regard  the  application  is 
utilized  as  a  tool  to  capture  the  relevant  physics  without  resolving  the  critical  is¬ 
sues  mentioned  above.  It  is  assumed  that  the  Baldwin-Lomax  model  will  perform 
adequately  in  the  applied  regime.  There  is  some  justification  for  this  assumption, 
as  the  pertinent  issue  is  the  determination  of  the  onset  of  unsteady  motion.  No 
statement  is  made  about  the  relevance  of  the  turbulence  model  after  a  dynamic  flow 
has  evolved.  That  is,  the  turbulence  model  does  not  necessarily  have  to  accurately 
model  the  large-scale,  unsteady  flow  structures  mentioned  above. 

The  turbulence  model  is  applied  in  three  separate  regions:  the  boundary  layer 
on  the  airfoil,  the  wake  proximal  to  the  airfoil  trailing  edge,  and  the  far-field  wake 
(see  Figure  22).  The  boundary  layer  of  the  airfoil  is  divided  into  inner  and  outer 
regions.  In  the  inner  region,  the  eddy  viscosity,  £,,  is  given  by  the  Prandtl-Van  Driest 
formula  [13] 

£<  =  p{KY D)2  M  ,  (64) 
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where  u>  is  the  vorticity,  Y  represents  the  normal  distance  from  the  airfoil  surface, 
K  =  0.40  is  the  von  Karman  constant,  and  the  Van  Driest  damping  factor  is  given 


by 


D  —  1  —  exp 


(65) 


The  subscript  s  denotes  conditions  at  the  surface  (wall). 

In  the  outer  region  of  the  boundary  layer,  the  eddy  viscosity  is  defined  by 


£o  —  pKclCcpYmaxFmaxFki  (66) 

where  Fmax  =  max(K \u\D)  ,  and  Ymax  is  the  value  of  Y  corresponding  to  Fmax.  The 
function  Fmax  is  maximized  over  the  Y  values  at  a  particular  station.  Values  for  the 
constants  are  taken  to  be  those  for  an  equilibrium  incompressible  boundary  layer 
[14],  for  which  Cc  =  1.2,  Ck  =  0.65.  The  intermitancy  factor,  F\,  is  given  by 

[l+5.5(Cfcy/rmax)6]'1,  (67) 

and  the  Clauser  constant,  K =  0.0168. 

The  turbulence  model  switches  from  the  inner  to  outer  formulation  at  the  first 
value  of  Y  for  which  e,  >  eQ.  Transition  from  laminar  to  turbulent  flow  is  specified 
by  user  input. 

The  turbulence  model  in  the  far-wake  is  modeled  by  defining  the  turbulent 
eddv  viscosity  there  as  [39] 

Y  At/2 

e*k  =  - Fk,  (68) 

*  max 

where 


Cwk  —  0.058, 
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Fma  x  =  max(V'|o;|), 

At;  =  («’ +  .»«£- («>  +  »»)!&. 

In  the  wake,  Ymax  is  measured  from  the  wake  centerline  as  determined  by  the 
location  of  minimum  velocity.  The  intermitancy  factor,  Fk,  is  again  obtained  from 
Eq.  (67).  The  value  of  the  constant  Cwk  was  chosen  to  match  the  theoretical  value 
for  an  incompressible  turbulent  wake  [40]. 

In  the  near- wake,  the  turbulent  eddy  viscosity  is  determined  by  algebraically 
smoothing  the  eddy  viscosity  profile  at  the  trailing  edge  to  the  fan-- wake  profile  [39]. 
The  distance  to  the  fair  wake  application  was  chosen  to  be  of  order  106,  where  6  is 
the  average  boundary  layer  thickness  at  the  trailing  edge  [39]. 


Boundar 
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y  Layer 

Near-Wake 

Far-Wake 

' 

< 

i 

** 

Xo 
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Figure  22.  Regions  of  Application  for  Turbulence  Model 
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4-2  Modification  of  the  Incompressible  Code 

In  order  to  compare  the  numericaliy  calculated  lift  with  the  lift  predicted  by 
Theodorsen’s  function,  it  is  necessary  to  forcibly  oscillate  the  airfoil  at  moderate 
rotational  rates  (reduced  frequencies,  k ,  greater  than  2).  The  incompressible  code 
utilized  in  Part  I  was  not  robust  enough  to  sustain  these  rotational  rates.  At  values 
as  low  as  k  =  2,  computational  instabilities  arose,  and  the  time  history  of  Ci  began 
to  display  asymmetric  behavior.  Difficulties  of  this  type  have  been  experienced  b^ 
other  investigators  [41,  42].  The  difficulties  arose  only  at  increased  rotational  rates. 
The  low  frequencies  associated  with  the  initial  onset  of  unsteady  motion  did  not 
destabilize  the  computations.  Hegna  [41]  and  Meh<a  [35]  have  also  applied  the  same 
procedure  with  success  when  the  rotational  rates  were  low. 

To  approach  the  higher  rotational  rates,  the  incompressible  code  was  modified 
to  employ  a  time-dependent  coordinate  transformation  [34]  of  the  form 

t  =  rj  =  T](x,y,t),  r  =  t,  (69) 

with  transformation  metrics  given  by 

£x  =  )  £il  =  ~XjjJ  1  it  =  ~Xrix  —  Vriy , 

T)x  =  -y(J~\  T)y  =  -X(j~\  T]t  =  -XtT)x  ~  yrVy 

The  grid  speeds  are  determined  by  xT  =  —  fly  and  yT  =  fix,  where  x  and  y  represent 
the  components  of  the  momen  arm  between  the  axis  of  rotation  and  the  point  in 
question  (see  Figure  23).  The  Jacobian  of  the  transformation  is  given  by 

j  —  Xtfr,  —  X^y^  7^  0. 
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Figure  23.  Definition  of  x,y  Coordinate  System 


Applying  this  transformation  results  in  an  identical  form  for  the  vorticity-transport 
equation  (cf.  Eq.  (29)) 


ut  +  uu>t  +  VU)V  =  -^V2u;, 
provided  the  convective  coefficients  are  modified  as 


(70) 


u  =  9T)J-1+(t,  i )  =  + 


In  a  similar  fashion,  the  calculation  of  pressure  presented  in  Part  I  (Eqs.  (34)-(35)) 
is  modified  to  produce 

(«*  +  +  m „)  +  -  ${«»))  +  j~\yr ,Pi  -  ViPv)  =  (71) 

(V<  +  fat  +  7?,V„)  +  J _1(^tJV€  -  »P*Vt,)  +  j_1(l{P^  -  X vPi)  =  (72) 

and  these  equations  are  manipulated  to  yield 

p<  =  i  _  +  ytvn)  -  +  iw)] 

-Xt(ut  +  (tUi  +  T}tu„)  -  yt(vt  +  fat  +  T)tv„),  (73) 
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(74) 


Pr,  =  J  -  <hv*)  +  Vt(XVUV  +  yvVV )  -  +  y*!^)] 

-x,(u<  +  £tUt  +  1 ituv)  -  yv(vt  +  ZtVt  +  T)t vv). 

The  quantities  fa  are  the  coefficients  associated  with  the  transformed  Laplace  oper¬ 
ator,  and  are  defined  in  Eq.  (28).  As  in  Part  I,  Eq.  (74)  is  integrated  along  a  line  of 
constant  £  to  obtain  the  pressure  at  the  leading  edge  of  the  airfoil,  while  Eq.  (73)  is 
integrated  along  the  upper  and  lower  airfoil  surfaces  to  obtain  the  pressures  there. 
No  further  simplification  of  Eq.  (73)  is  assumed  along  the  airfoil  surface. 

A  fundamental  difference  in  the  present  approach,  as  compared  with  the  method 
presented  Part  I,  is  the  application  of  the  boundary  conditions.  In  Part  I,  the  equa¬ 
tions  were  expressed  in  a  non-inertia!  frame,  and  the  inclusion  of  the  apparent  body 
forces  (centrifugal,  tangential,  and  Coriolis)  accounted  for  the  non-inertial  character 
of  the  reference  frame.  For  such  an  approach,  the  velocities  on  the  airfoil  surface  are 
zero  and  the  freestream  velocities  vary  as  the  airfoil  oscillates.  Conversely,  when  ap¬ 
plying  the  time-dependent  coordinate  transformation,  the  freestream  velocities  are 
fixed  and  the  surface  velocities  are  equal  to  the  actual  velocity  at  a  point  on  the 
moving  surface.  This,  of  course,  simplifies  the  application  of  the  freestream  veloci¬ 
ties,  but  complicates  those  applied  on  the  airfoil  surface.  The  non-zero  velocities  on 
the  surface  necessitate  a  change  in  the  definition  of  the  operator  L3  (cf.  Eq.  (57)) 
and  the  expression  for  the  discrete  boundary  condition  for  vorticity  on  the  airfoil 
surface 

L3Vi  +  +  S  =  0,  (75) 

where  S  represents  a  known  source  term  arising  from  the  prescribed  velocities  on 
the  airfoil  surface.  In  delta  form,  Eq.  (75)  becomes 

£,3 An 0  +  Anu>  =  -  L3tl>n  —  ion  —  L3yn+1  -  Sn+1 ,  (76) 
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where 


y  =  y  cos  a  —  xsina,  (77) 

and  the  components  of  5  are  given  by 

S,n+‘  =  j-2  [-2(*  +  *)a"+1  +  *<."«  +  *<+>  + 

+§«?,'  -  §«i'  +  ^cT+'  +  «+I] ,  (78) 

and  to  first-order  accuracy  in  time,  which  is  consistent  with  the  overall  time-accuracy 
of  the  scheme, 

a?+1  =  ( uy  -  ux)”+1  =  (tty  -  t>x)"  -  At  [fi(yy  +  xx)  -  S l2{xy  -  yx)]", 

6"+1  =  (2ux„  -  2ttyn)“+1  =  (2t>x„  -  2uy„)”  +  2A tfi(xx„  +  yy„)", 
c"+1  =  («y«  “  *>*€)"+1  =  -  vx()?  -  A ttl{yyi  +  xx()", 

<*?+1  =  («y^  “  ux,)7+1  =  ( uy„  -  yx,)^  -  A ffl(yy„  +  xx„)?.  (79) 

The  superscript  n  refers  to  the  temporal  index  while  the  subscript  i  refers  to  the  node 
index  on  the  airfoil  surface.  The  terms  x  and  y  represent  the  respective  moment  arms 
from  a  point  on  the  airfoil  surface  to  the  center  of  rotation.  The  terms  SI  and  SI  are 
known  from  the  prescribed  nature  of  the  rotation. 

In  addition,  the  streamfunction  on  the  surface  is  now  nonzero.  Thus  Eq.  (45) 
is  modified  to 


y.  =  'P.  +  y,  =  -^(yy  +  xx)  =  -Sir, 


(80) 


f  =  (yy  +  xx), 


leading  to 


+  y»  +  Sir  =  0. 


(81) 
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In  delta  form,  the  surface  boundary  condition  on  streamfunction  is  written  as 


AntP,  =  -(Vf  +  y,n)  -  A ny,  -  nn+1  P+1 ,  (82) 

with  first-order  accurate  approximations 

=  yn+1  cos  a  —  xn+1  sin  a  —  yn  cos  a  +  xn  sin  a 
=  (yn+1  —  yn)  cos  a  —  (xn+1  —  xn)  sin  a, 

Any,  =  At(yrcosa  —  xTsina),  (83) 

and 

1T+1  =  +  A  thn,  (84) 

fn+1  =  (yy  +  xx)n+1, 

fn+1  =  (yn  +  yTA t)(yn  -f  yTA<)  -f  ( xn  +  xTA/)(z"  +  xTA t).  (85) 

A  comparative  Fourier  (von  Neumann)  stability  analysis  between  the  approach 
indicated  in  Part  I  and  the  time-dependent  metric  approach  described  in  this  section 
is  presented  in  Appendix  B.  The  stability  analysis  is  applied  to  the  linearized,  viscous 
Burger’s  equation  modeling  the  vorticity-transport  equation. 

4-3  The  Genesis  of  Theodorsen’s  Function 

The  unsteady  aerodynamic  forces  on  an  oscillating  airfoil  in  potential  flow  were 
determined  by  Theodorsen  [15]  in  1935.  The  development  is  based  on  potential  flow 
theory  and  an  assumed  distribution  of  circulation  in  the  wake.  The  magnitude  of 
the  circulation  is  established  by  imposing  the  Kutta  condition  at  the  airfoil  trailing 
edge.  The  solution  is  resolved  into  Bessel  functions  of  the  first  and  second  kind, 
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and  is  expressed  in  terms  of  the  reduced  frequency  k.  Small  perturbations  about  an 
equilibrium  solution  are  applied  and  the  resultant  motion  is  assumed  harmonic.  The 
development  of  the  velocity  potentials  proceeds  by  dividing  them  into  two  groups: 
those  due  to  the  noncirculatory  flow,  and  those  due  to  circulatory  flow. 

Pursuing  the  noncirculatory  potentials  first,  the  airfoil  is  represented  by  a  unit 
circle  (see  Figure  24).  The  potential  of  a  source  of  strength  e  at  the  origin  is  given 
by  [15] 

<P  =  ~  log(z2  +  y2).  (86) 

For  a  source  c  at  (xi,yi)  on  the  unit  circle,  the  potential  is 


'p  =  ir  logK*  ”  **)’  +  l-y~  vi)2l- 


(87) 


Placing  a  double  source  2c  at  (ii,yi)  and  a  double  negative  source  at  (xj,  — yi),  the 
potential  becomes 

2 *  8(*-xi)2  +  (y  +  yi)2 
The  function  <p  on  the  unit  circle  provides,  upon  integration,  the  surface  potential 
of  line  pq  (see  Figure  24).  To  determine  the  potential  due  to  a,  substitute  for  the 
source  strength 

e  =  -Uab,  (89) 

where  U  is  the  freestream  velocity,  and  b  is  the  half-chord.  Taking  y  =  y/l  —  x2  on 
the  unit  circle  and  integrating  yields  the  potential  due  to  angle  of  attack, 


¥>a  = 


—Uab 

2n 


(x  —  xi)2  +  (y  —  yi)2  , 
(x  -  Xi)2  +  (y  +  yi)2  Xl 


(90) 


or 


V?a  =  Uaby/l  —  x2. 


(91) 
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The  potential  due  to  an  airfoil  in  downward  motion  with  a  velocity  h  (positive 
downward)  is  obtained  by  substituting  h/U  for  a  in  Eq.  (91). 

<fik  =  hbV  1  -  x2.  (92) 

The  potential  for  an  airfoil  rotating  about  a  point  located  a  distance  a  (positive  to 
the  right)  from  the  midchord  point  is  determined  by  considering  a  superposition  of 
potentials.  The  first  being  the  potential  associated  with  rotation  about  the  leading 
edge,  and  the  second  associated  with  the  vertical  motion  with  velocity  — d(l  +  a)b. 
The  latter  results  in  a  modification  to  Eq.  (92),  resulting  in 

(fi-k  =  — d(l  +  a)b2y/ 1  —  x2.  (93) 


The  potential  due  to  rotation  about  the  leading  edge  with  angular  velocity  a  is 
described  by  setting  e  =  — (xi  +  l)d62  in  Eq.  (88)  and  integrating 


<PalB  — 


(x  -  xi)2  +  (y  -  yi)2 
(x-x!)2  +  (y  +  yi)2 


(x,  +  l)dx,, 


(94) 


or 

£[, 2 

<PalE  =  ~2~(X  +  2Wl  ~  x2‘  (95) 

As  indicated,  the  potential  for  an  airfoil  rotating  with  angular  velocity  a  about  a 
prescribed  axis  of  rotation  is  the  sum  of  Eqs.  (93)  and  (95) 


=  <Pa  le  +  (96) 

ipo,  =  ab2(\x  -  a)Vl  —  x2.  (97) 

Employing  the  extended  Bernoulli  equation,  the  local  pressure  is,  to  within  a  spatially 
invariant  function  of  time 

p  =  -|(v'I+2w)’  (98) 


58 


Figure  24.  Conformal  Map  of  Wing  Profile  with  Circulation  Element  in  Wake 


where  V  is  the  local  velocity  and  <p  the  velocity  potential  at  a  point.  Substituting 
V  =  U+jg,  the  pressure  difference  between  the  upper  and  lower  surfaces  at  a  point 
x  is  given  by 


Substituting  the  individual  velocity  potentials  (Eqs.  (91)  and  (97))  into  Eq.  (99) 
and  integrating  produces  the  net  force  on  the  airfoil  for  the  noncirculatory  potentials 


Po  = 


—pb2[Uira  —  afore*]. 


(100) 


In  a  like  manner,  the  circulatory  velocity  potentials  due  to  an  assumed  distri¬ 
bution  of  circulation  of  strength  A(x)  extending  from  the  trailing  edge  to  i  =  +oo 
(see  Figure  24)  are  derived.  The  potential  at  a  point  (x^yj)  on  the  unit  circle  due 
to  a  vortex  element  pair  of  strength  A T  is  given  by  Theodorsen  [15] 


¥>r(xj,yi)  = 


AT 

27T 


arctan 


-  “cUn  (iTTv)]  ’  (101) 
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Introducing  the  transformation 


Xo  +  X;1  =  2x0,  (102) 

or  on  the  x  axis 

Xo  =  Xo  +  y/xl  -  1,  (103) 

with  the  further  substitutions  (on  the  unit  circle) 

Xi  =  x,  (104) 

yi  =  Vl  -x2,  (105) 


provides 


(106) 


Equation  (106)  gives  the  clockwise  circulation  about  the  airfoil  due  to  the  vortex 
element  —AT  at  xa.  The  extended  Bernoulli  equation  (Eq.  (99))  is  again  applied. 
However,  the  vortex  element  —  Ar  is  regarded  as  converting  to  the  right  relative  to 
the  airfoil  with  velocity  U 


dy  _  dy  dx0  _  dtp 
dt  dx0  dt  dx0 


Hence  Bernoulli’s  equation  is  written 


p  =  —2  pU 


Differentiation  of  Eq.  (106)  provides 

dp  _  AY  /eq-1 

dx  2tt  (x0  —  x)y/l  —  i2’ 


(107) 


(108) 


(109) 
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and 


8<p 

dxa 


Ar  y/T- 


2*  (x0  -  x)yjxl  -  1 
Integration  provides  the  force  on  the  airfoil  due  to  the  element  — Af 


(110) 


APr  =  -pUbAT 


(111) 


(112) 


Further  integration,  with  Ar  =  A dx0  along  the  wake,  provides  the  total  circulatory 
force  on  the  airfoil 

Pr  =  -pUb  r  —f==S.dx0.  (113) 


The  magnitude  of  the  circulation  is  determined  by  imposing  the  Kutta  condi¬ 
tion,  requiring  that  the  velocity  components  are  finite  at  the  trailing  edge 


+  Vi  +  Vr)  /  oo. 


(114) 


Introducing  the  velocity  potentials  (Eqs.  (91), (97),  and  (106)),  Eq.  (114)  is  further 
expressed  by  the  relationship 


Q  =  ~  r  Adx0  =  Ua  +  b(\-  a)a.  (115) 

Zir  J\  v  x0  —  1 

The  distribution  of  circulation  in  the  wake  is  assumed  harmonic,  with  the  substitu¬ 
tion 


A  =  A0exp{i[k(3/b-  xa)  +  $},  (116) 

where  s  =  Ut  is  the  distance  from  the  first  vortex  element  to  the  airfoil,  i  =  \/— T, 
and  k  is  the  reduced  frequency  representing  the  wavelength.  The  circulatory  forces 
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Figure  25.  Real  and  Imaginary  Components  of  the  Theodorsen  Function 
on  the  airfoil  then  become 


Pr  =  —2pUbnQC(k), 


(117) 


where  the  complex  function  C(k )  is  the  Theodorsen  function,  given  by 


j:*? 

C(k)  =  J  1  v  ° 


-iJeic 


dxa 


/OO 


t~'kx°dxn 


(118) 


The  real  and  imaginary  parts  for  the  complex  function  C(k)  =  F(k)  +  iG(k) 
are  shown  in  Figure  25. 
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The  net  force  on  the  airfoil  is  then  obtained  by  summing  the  forces  due  to  the 
noncirculatory  potentials  with  those  of  the  circulatory  potentials 

P  =  Po  +  Pr, 

=  —  pvb2(Ua  —  aba)  —  2vpUbC{k)  |t/a  +  6(|  —  a)dj .  (119) 

The  expression  for  P  is  then  nondimensionalized.  Along  with  a  small  angle  as¬ 
sumption  and  the  substitutions  0  =  —a,  0  =  — ( c/U)a ,  and  0  =  — ( c/U)2a ,  the  lift 
coefficient  can  be  expressed  as 

Ck  =  \*0  +  \air§  —  C(k)  +  (|  —  a)x0j  .  (120) 

where  the  time  derivatives  of  0  are  now  taken  with  respect  to  nondimensional  time, 
i.e.,  aerodynamic  time  units  nondimensionalized  by  c/U.  This  is  the  form  used  for 
comparison  with  the  numerical  results.  It  should  be  noted  that  for  the  steady  case 

M-0,  (121) 

with 

F(k)  — *  1  and  G(k)  —*  0,  (122) 

so  that 

Ck  =  -2n0  =  2ira.  (123) 

as  predicted  by  thin-airfoil  theory  [43]. 

4-4  Beam  and  Warming  Code 

The  Beam  and  Warming  code  of  Visbal  [27]  applied  to  the  two-dimensional, 
compressible  Navier-Stokes  equations  is  utilized  for  comparison  with  the  modified 
incompressible  code  and  the  results  predicted  by  Theodorsen’s  function.  The  Beam 


63 


and  Warming  algorithm  [44]  is  an  implicit  approximate-factorization  scheme,  and  is 
implemented  by  Visbal  for  a  moving  O-grid  configuration.  The  same  time-dependent 
coordinate  transformation  described  earlier  is  utilized 


£  =  v  =  v(x,yJ),  t  =  t. 


The  resulting  equations,  in  strong-conservation  law  form  [39],  are 

dU  dEi  dE2  _ 

at  +  as,  +  art 

mvA)  av2{u,v„)  aw,(u,v()  dw2(u,u„ > 

d(  +  d(  +  dr,  +  ft) 


(124) 


where 


pU 

pvU  +  (xp 
pvU  +  ZyP 
(p  +  pe)l4  -  hp 


E2 


pV 

puV  +  T)xp 
PVV  -I-  TjyP 

(p  +  pe)V  -  t itP 


(125) 


(126) 


(127) 
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0 

+  bjVt 

biting  +  ^(vu^  -f-  uv^)  +  1>3VV(  +  b4T( 

0 

Cl  Ur,  +  C2VV 
C3U„  +  c4v„ 

CiUUr,  +  C2UVr,  +  C3UU„  +  CjWt?,,  -(-  C5T„ 

0 

Cl  U£  +  C3V( 
c2u(  +  C4t>£ 

cxuu^  +  C2VU£  +  c3urf  +  c4i>v£  +  c5T£ 

0 

+  d2v„ 
d2uv  +  d3vv 

diUUr,  +  d2(vuv  +  uvv)  +  d3vvv  +  d4Tv 

where  U  and  V  denote  the  contravariant  velocities 

W  =  ixu  +  iyv  +  &,  (132) 

V  =  7?*u  +  T)yV  +  T}t.  (133) 

and  6, ,  c,,  and  d,  represent  the  viscous  coefficients  [39],  Closure  for  the  equations 
is  enforced  by  Sutherland’s  viscosity  formula,  a  constant  Prandtl  number,  and  the 
perfect-gas  law. 

The  implicit  implementation  is  applied  to  the  strong  conservation  form  of  the 
Navier-Stokes  equations.  In  delta  form,  with  first-order  Euler  time  differencing,  the 
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scheme  is  written  as  [14]: 


-At  |^(£i  -  V,  -  V2)n  +  |-(£2  -Wt-  W2)"  J . 


(134) 


The  scheme  provides  the  correction  vector,  A U,  to  the  current  solution,  U " 


_  (jn  +  Af) ^ 


(135) 


where  n  represents  the  temporal  index  such  that  t  =  (n  —  1)A/  and  .4,  F,  M,  and 
N  represent  the  Jacobian  matrices  [39]. 


dEi 

dU' 
dV , 

dU,:' 


(136) 

(137) 


To  maintain  numerical  stability  and  provide  smooth  solutions,  explicit  fourth- 
order  damping  is  added  to  the  right-hand  side  of  Eq.  (134) 


D(  =  - f'AtJ~j6*U £, 


(138) 


and  implicit  fourth-order  damping  is  applied  to  the  left  side  of  Eq.  (107) 


A,  =  -/,A  tJ-^U^ 


(139) 


66 


where  6  is  a  second-order  accurate,  central-difference  operator.  The  recommended 
values  of  the  damping  coefficients  are  [45] 

fe  =  0(0.01),  /,  =  2/e,  (140) 

The  values  employed  for  this  application  are 

/<=  =  0.02,  /,  =  0.04.  (141) 
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V.  Results  and  Conclusions:  Part  II 


5. 1  Overview 

The  focus  of  Chapter  5  is  delineated  into  two  major  areas.  The  first  area  of 
exploration  is  the  numerical  prediction  of  flutter  onset  in  the  turbulent  regime.  As 
indicated  in  Part  I  (cf.  Figure  1),  there  is  a  fundamental  difference  between  the  lam¬ 
inar  and  turbulent  regimes,  and  this  is  reflected  in  the  numerical  results.  The  onset 
of  flutter  in  the  laminar,  low-speed  regime  occurs  at  a  significantly  lower  Reynolds 
number  than  the  experimental  results  indicate  [11,  15].  The  original  rationale  for  in¬ 
corporating  a  turbulence  model  was  to  address  this  issue.  The  incompressible  code 
is  well  validated  in  the  low-speed,  laminar  regime  (see  Appendix  A)  for  the  case 
of  the  circular  cylinder.  However,  the  experimental  efforts  for  the  airfoil  case  have 
been  exclusively  concentrated  in  the  higher  Reynolds  number  range  [Re  >  5  x  105) 
[11, 15,  12].  Therefore,  the  Baldwin-Lomax  turbulence  model  is  incorporated  to  cap¬ 
ture  the  experimental  data.  Previous  computational  efforts  in  this  regime  include 
the  work  of  Strganac,  Mook  and  Mitchum  [21],  who  numerically  investigated  sub¬ 
sonic  flutter  for  a  finite  wing  using  a  potential  flow  model.  As  previously  discussed, 
Guruswamy  [16]  completed  a  similar  effort  using  the  full  Navier-Stokes  equations. 
Robinson  et  ad.  [20]  have  applied  the  Euler  equations  to  predict  the  aieroelastic 
behavior  of  a  wing  with  a  deforming  mesh. 

The  second  area  examined  in  Chapter  5  is  the  correlation  between  the  numer¬ 
ically  predicted  aerodynamic  tramsfer  function  and  that  obtained  by  Theodorsen’s 
potential  flow  development.  The  effort  is  undertaken  with  the  goad  of  establishing  a 
basis  for  the  numericad  prediction  of  flutter  onset.  This  is  accomplished  by  compar¬ 
ing  the  theoreticad  lift  coefficient  predicted  by  Theodorsen’s  function  to  the  numer- 
icadly  predicted  lift  coefficient.  A  NACA  0012  aurfoil  is  oscillated  with  a  prescribed 
frequency  and  amplitude  to  generate  the  unsteady  lift  coefficient.  The  modified  in¬ 
compressible  code  and  the  compressible  code  of  Visbad  are  applied  in  this  capaicity. 
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Chaing  and  Fleeter  [46]  have  undertaken  a  similar  computational  effort,  computing 
the  unsteady  pressure  distribution  over  an  airfoil  and  comparing  the  results  to  those 
generated  by  Theodorsen’s  method.  This  was  accomplished  by  applying  a  locally 
analytical  method  to  the  solution  of  the  potential  flow  equations. 

5.2  Flutter  Onset:  Comparison  with  Theory  and  Experiment 

In  this  section  results  a re  presented  in  which  the  numerically  predicted  flutter 
onset  velocity  and  those  established  by  experiment  and  predicted  by  theory  are 
compared.  The  numerical  results  axe  obtained  using  Grid  1,  with  a  NACA  0015 
profile  substituted,  as  per  the  experiments  of  Yang  and  Zhao  [11].  The  incompressible 
code  used  in  Part  I  is  used  along  with  the  modified  Baldwin-Lomax  turbulence  model. 
To  allow  operation  of  the  code  at  the  higher  Reynolds  number  range,  77-upwinding 
of  identical  form  to  the  ^-upwinding  described  in  Part  I  is  applied.  The  theoretical 
and  additional  experimental  results  are  due  to  Theodorsen  [15]. 

The  numerical  procedure  employed  is  similar  to  the  method  described  in  Part  I 
for  the  laminar  case.  A  steady  solution  is  obtained  at  relatively  low  Reynolds  number 
(Re  =  2  x  105).  The  Reynolds  number  is  then  incremented  and  a  new  solution 
is  obtained.  The  increase  in  Reynolds  number  introduces  perturbations  into  the 
numerical  flowfield,  which  are  either  damped  or  amplified  to  an  unsteady  solution. 
The  Reynolds  number  increase  was  on  the  order  of  5  x  104,  however,  smaller  incre¬ 
ments  (~  104)  were  applied  when  the  aerodynamic  damping  was  reduced,  i.e.,  as  the 
critical  point  was  approached. 

The  results  are  shown  in  Figure  26.  The  general  trend  of  the  numerical  data  is 
correct.  However,  the  results  are  seen  to  be  a  function  of  the  natural  frequencies,  u>0 
and  u>h,  rather  than  the  ratio.  This  is  in  contrast  to  the  theoretically  predicted  result 
[15].  The  reason  for  this  effect  is  the  subcritical  nature  of  the  bifurcation.  At  lower 
values  of  u>a  and  wj,,  a  more  dynamic  augmented  system  is  provided;  perturbations 
are  amplified  more  readily  in  this  case.  The  presence  of  perturbations  of  increased 
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magnitude  then  allow  the  possiblity  of  approaching  a  limit-cycle  attractor  rather 
than  an  equilibrium  attractor.  As  a  further  consideration,  consider  a  pitch-axis 
case  where  the  value  of  u>Q  is  fixed,  but  Ka,Ia  <C  1-  In  this  case  it  is  obvious 
that  the  airfoil  will  be  extremely  responsive  to  perturbation.  Additionally,  Yang 
and  Zhao  [11]  have  observed  a  plunge-dominated,  unsteady  mode  that  could  only 
be  excited  by  providing  an  initial  displacement  in  plunge.  This  result  represents  a 
differential  response  to  perturbation  level  and  thus  offers  experimental  evidence  of 
subcritical  bifurcation.  The  results  shown  in  Figure  26  are  explained  on  this  basis, 
i.e.,  subcritical  bifiurcation  and  a  consequent  differential  response  to  perturbation 
level.  For  instance,  as  u>a  is  decreased  (for  a  fixed  airfoil  mass  and  momerat  of 
inertia),  the  torsional  spring  constant  Ka  is  reduced.  The  weaker  torsional  spring 
allows  a  more  dynamic  response  to  aerodynamic  perturbation,  i.e.,  the  aerodymamic 
perturbations  introduced  by  the  Re- continuation  method  explained  above. 

There  are  additional  sources  of  error  in  the  numerical  computations  which 
should  be  addressed.  Firstly,  the  Baldwin-Lomax  turbulence  model  has  been  criti¬ 
cized  in  unsteady  flow  applications  [14,  47],  and  in  the  ability  to  predict  self-induced 
oscillations  [16].  Secondly,  compressibility  effects  remain  unresolved  and  are  likely 
a  factor  in  this  regime.  Thirdly,  when  employing  the  turbulence  model  at  higher 
Reynolds  numbers,  timestep  requirements  became  too  restrictive,  and  therefore,  ex¬ 
amining  results  with  a  finer  grid  was  computationally  prohibitive.  And  lastly,  in  the 
course  of  subsequent  computations,  it  became  apparent  that  the  concurrent  appli¬ 
cation  of  77-  and  £-upwinding,  when  combined  with  the  non-conservative  form  of  the 
equations,  impacted  the  phase  response  of  the  aerodynamic  output.  Nonetheless,  the 
results  do  indicate  clearly  the  demarcation  between  these  results  and  those  presented 
in  Part  I,  and  the  role  of  turbulence  in  this  demarcation  (See  Figure  1). 
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Figure  26.  Flutter  Onset  Velocity  as  a  Function  of  u>h/va. 
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5.3  Comparison  with  Inviscid  Aerodynamic  Transfer  Function 

In  this  section,  a  basis  is  established  for  the  application  of  the  considered 
numerical  algorithms  to  the  prediction  of  flutter  onset  for  an  airfoil.  The  theoretical 
development  of  the  prediction  of  flutter  was  pioneered  by  Theodorsen  [15],  using 
potential  flow  theory  in  conjuntion  with  an  assumed  distribution  of  circulation  in 
the  wake  and  the  application  of  the  Kutta  condition  at  the  trailing  edge.  In  the 
theoretical  approach,  the  derivation  of  Theodorsen’s  function  models  the  inviscid 
aerodynamic  transfer  function.  It  provides  the  aerodynamic  response  to  an  input  of 
prescribed  amplitude  and  frequency. 

The  general  intent  of  the  numerical  efforts  in  this  regard  is  to  compare  the 
aerodynamic  transfer  function  described  by  Theodorsen  to  that  predicted  numeri¬ 
cally.  This  is  done  for  two  codes:  the  modified  form  of  the  incompressible  code  and 
the  compressible  code  of  Visbal  [27],  modeling  the  full  Navier-Stokes  equations  (see 
Chapter  4).  The  approach  is  to  apply  a  forcing  function,  angle  of  attack,  a(t),  of 
known  amplitude  and  frequency  and  assess  whether  the  numerical  results  correspond 
to  the  theoretical  results  of  Theodorsen,  across  a  broad  range  of  frequencies.  Since 
Theodorsen’s  approach  is  derived  assuming  small  perturbations  about  an  equilibrium 
state,  the  amplitude  of  the  angular  oscillation  is  fixed  at  0.01  radians.  The  airfoil 
is  oscillated  at  reduced  frequencies,  k ,  ranging  from  0.05  to  6,  with  the  Reynolds 
number  varied  between  Re  =  100  and  Re  =  3  x  106.  All  cases  axe  initiated  from  an 
equilibrium  solution  at  q  =  0  for  a  given  Reynolds  number. 

Results  are  compared  in  Figures  27-32  for  the  compressible  code  and  Fig¬ 
ures  33-36  for  the  incompressible  code.  The  lift  coefficient  predicted  by  Theodorsen’s 
function  is  denoted  C*,  (cf.  Eq.  (120))  and  is  represented  by  long-dashed  lines  in  Fig¬ 
ures  27-36.  Ci  refers  to  the  numerically  computed  lift  coefficient,  and  is  represented 
by  solid  lines.  The  angle  of  attack,  a,  is  indicated  by  short-dashed  lines. 

The  agreement  with  Theodorsen’s  function  is  seen  to  be  a  function  of  the 
reduced  frequency  k ,  but  in  general  the  results  are  good.  An  exact  correspondence 
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was  not  anticipated;  differences  in  the  approaches  remain,  primarily  consisting  of: 

•  Reynolds  number  effects 

•  airfoil  thickness  effects 

•  validity  of  the  assumed  distribution  of  circulation 

•  grid-related  effects 

•  artificial  dissipation  effects 

•  compressibility  effects 

The  theoretical  result  is  derived  assuming  an  inviscid  fluid,  corresponding  to  an 
infinite  Reynolds  number.  The  numerical  computations  are  constrained  to  finite 
Reynolds  numbers,  and  therefore  some  differences  due  to  Reynolds  number  effects 
may  remain.  This  is  particularly  true  for  the  incompressible  code,  since  it  is  con¬ 
strained  to  a  lower  Reynolds  number  range  than  the  compressible  code.  The  dis¬ 
tribution  of  circulation  assumed  in  the  theoretical  development  (cf.  Eq.  (116))  in 
conjunction  with  the  application  of  the  Kutta  condition  (cf.  Eq.  (115)),  allows  a 
framework  for  incorporating  viscous  effects  into  the  inviscid  analysis.  The  manifesta¬ 
tion  of  vorticity  in  the  wake  and  the  satisfaction  of  the  Kutta  condition  for  the  viscous 
codes  are  accomplished  by  a  completely  different  mechanism,  i.e.,  the  direct  solution 
of  the  field  equations  with  am  associated  set  of  boundary  conditions.  Grid-related 
effects  are  important  because  the  moving  grid  establishes  time-dependent  numerical 
errors  which  must  be  damped.  The  ability  to  damp  these  errors  is  dependent  upon 
the  grid  resolution.  Errors  of  this  type  occur  in  addition  to  the  discretization  or 
truncation  errors  associated  with  numerical  computations  performed  on  a  fixed  grid. 
Artificiad  dissipation  can  further  influence  the  numerical  results,  for  instance,  it  was 
mentioned  eairlier  that  concurrent  application  of  £-  and  r/-upwinding  introduces  phase 
errors  into  the  numerical  calculations.  And  lastly,  for  the  case  of  the  compressible 
code,  compressibility  effects  are  a  factor.  Of  the  above  items,  only  Reynolds  number 
effects  die  addressed  in  a  systematic  fashion. 
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The  most  meaningful  trend  is  the  convergence  of  the  numerical  results  to  those 
of  theory  as  the  Reynolds  number  increases,  i.e.,  as  an  inviscid  solution  is  approached. 
This  trend  is  evident  at  every  frequency  considered,  except  for  one  case  for  the  in¬ 
compressible  code,  which  will  be  addressed  later.  The  results  for  the  compressible 
code  are  discussed  first.  More  runs  were  performed  with  this  code  because  the  re¬ 
quired  modifications  were  much  simpler,  and  data  was  collected  during  the  develop¬ 
ment  of  the  modified  incompressible  code.  All  displayed  results  for  the  compressible 
code  were  performed  using  an  O-grid  with  grid  characteristics  equivalent  to  Grid  1. 
The  freestream  Mach  number  was  set  at  0.2,  with  fourth-order  explicit  and  implicit 
damping  levels  fixed  at  fe  =  0.02  and  /,  =  0.04.  No  second-order  damping  was 
employed.  Previous  numerical  studies  [48,  49]  indicate  that  the  timestep  required 
for  an  accurate  temporal  resolution  may  be  less  than  that  required  for  numerical 
stability.  To  ensure  accuracy,  the  timestep,  A t,  is  adjusted  to  enforce  a  minimum  of 
1000  timesteps  per  period  of  oscillation,  as  recommended  by  Stanek  and  Visbal  [49]. 

The  first  result  shown  is  for  k  =  0.05,  see  Figure  27.  As  the  Reynolds  number 
is  increased  from  100  to  3  x  10®,  a  clear  convergence  to  theory  is  apparent.  The  trend 
occurs  for  each  increment  in  Reynolds  number;  the  amplitude  appears  almost  exact 
for  the  Re  =  3  x  106  case.  A  slight  phase  error  between  C;  and  C*  is  not  rectified 
by  the  increased  Reynolds  number.  The  k  =  0.2  case  shows  an  identical  trend;  the 
amplitude  is  nearly  exact,  but  the  phase  lag  is  slightly  exacerbated  at  this  frequency. 
At  k  =  1,  an  interesting  transition  appears,  as  the  error  appears  to  express  itself  as 
an  amplitude  rather  than  a  phase  error.  The  amplitude  shows  a  slight  error  while 
the  phase  lag  is  negligible.  An  increase  to  k  =  2  appears  to  resolve  both  the  phase 
and  amplitude  errors,  as  the  solution  at  Re  =  3  x  106  displays  almost  no  amplitude 
or  phase  error.  The  results  for  k  =  4  mimics  the  results  for  k  =  0.05  in  that  the 
amplitude  error  now  appears  negligible  but  the  phase  error  is  reintroduced.  Finally, 
the  results  for  k  =  6  appear  to  follow  the  same  trend  as  the  k  =  0.2  case,  i.e., 
the  phase  lag  is  exacerbated  compared  with  the  prior  result.  The  results  indicate  a 
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cyclic  trend  in  the  phase  and  amplitude  errors.  The  graphical  presentation  of  this 
trend  is  shown  in  Figure  38,  which  summarizes  the  phase  and  amplitude  errors  as  a 
function  of  k,  where  the  phase  error  is  indicated  as  a  percentage  of  the  period  and  the 
amplitude  error  is  based  on  peak  values.  It  is  expected  that  beyond  k  =  6,  numerical 
errors  would  begin  to  predominate,  and  that  poor  numerical  results  would  eventually 
ensue.  See  the  stability  analysis  presented  in  Appendix  B  for  am  indication  of  this 
behavior. 

The  results  for  the  incompressible  code  are  displayed  in  Figures  33-36.  Several 
difficulties  arose  in  the  application  of  the  incompressible  code.  The  primary  difficulty, 
already  mentioned,  is  the  problem  associated  with  concurrent  application  of  £-  and  77- 
upwinding  in  conjunction  with  the  non-conservative  form  of  the  equations.  It  became 
apparent  that  the  application  of  77-upwinding  produced  unwanted  phase  errors  in 
the  solution.  However,  77-upwinding  is  required  to  approach  the  higher  Reynolds 
numbers.  Therefore,  77-upwinding  was  not  applied  to  assess  the  correspondence  with 
Theodorsen’s  function,  This,  in  turn,  imposed  a  restriction  on  the  Reynolds  number 
range,  and  a  more  stringent  timestep  requirement.  For  this  reason,  the  maximum 
Reynolds  number  considered  in  this  case  was  2000.  However,  this  proved  satisfactory 
to  demonstrate  the  trends  previously  observed  with  the  compressible  code.  It  should 
be  noted  that  for  Re  >  4000,  the  numerical  results  were  degraded,  as  a  function  of 
k ,  for  the  reasons  expressed  above. 

Figure  33  shows  the  progression  in  Reynolds  number  for  k  =  0.05.  A  clear 
convergence  to  the  theoretical  result  is  indicated.  However,  at  Re  =  2000  there 
is  still  a  small  amplitude  error.  Figure  34  displays  the  results  for  the  k  =  0.2 
case.  Again,  a  clear  convergence  is  noted,  with  a  very  close  agreement  in  both 
phase  and  amplitude  for  Re  =  2000.  The  results  for  k  =  2  (Figure  35)  also  show 
close  agreement  at  Re  =  2000.  The  limitations  of  the  incompressible  code  becomes 
apparent  at  higher  frequencies,  as  the  results  for  k  =  4,  shown  in  Figure  36,  indicate. 
The  results  appear  to  be  converging  as  Re  increases  from  100  to  500,  however,  when 
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Re  is  increased  to  1000,  the  solution  is  degraded  rather  than  improved.  This  result  is 
attributed  to  the  lack  of  dissipation  in  the  r?  direction.  This  argument  is  supported 
by  the  apparent  convergence  for  the  lower  Reynolds  numbers,  where  the  natural 
dissipation  is  presumably  sufficient.  However,  as  the  Reynolds  number  continues 
to  increase,  the  natural  dissipation  becomes  insufficient.  Additionally,  this  problem 
predominates  only  at  higher  frequencies,  when  the  oscillations  in  velocity  normal 
to  the  surface  assume  a  larger  magnitude,  exacerbating  the  difficulties  due  to  the 
lack  of  77-upwinding.  The  trend  continues  to  develop  as  the  solution  at  Re  =  2000 
exhibits  increased  amplitude,  with  slight  variation  in  peak  amplitudes  noted. 

The  convergence  to  the  theoretical  result  is  quantified  by  defining  a  norm 

IIACil  =  ^(C,  -  C*)2/*,  (H2) 

where  N  is  the  number  of  discrete  points  over  which  the  summation  is  applied.  The 
summation  is  applied  over  1  periods,  where  /  is  an  integer.  The  beginning  of  a 
period  was  taken  to  be  when  Ck  =  0,  which  generally  corresponded  to  the  least 
error.  The  variation  of  ||AC/||  with  Re  for  k  =  0.2  is  displayed  in  Figure  39  for  the 
incompressible  code. 

A  parametric  investigation  is  presented  in  Table  3.  The  results  were  computed 
using  the  compressible  code,  at  k  =  1.  Case  1  represents  the  baseline  result  presented 
in  Figure  29d.  In  cases  2  through  4  the  artificial  dissipation  levels,  fe  and  were 
reduced  by  a  factor  of  two,  successively.  In  cases  2  and  3  there  is  a  slight  increase  in 
||AC/||  ,  indicating  that  the  initial  dissipation  levels  were  required  to  provide  suffi¬ 
cient  damping  for  the  case  of  the  oscillating  airfoil.  The  dissipation  levels  employed 
in  case  4  were  insufficient  and  resulted  in  divergence.  Case  5  explores  the  influence 
of  the  turbulence  model.  With  the  turbulence  model  active,  there  is  no  discern- 
able  change  in  the  computed  results,  and  ||AC/||  remains  invariant  to  five  significant 
figures.  Cases  6  through  11  explore  the  impact  of  compressibility  by  incrementing 
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the  Mach  number  from  0.2  to  0.5.  Initially,  (case  6)  there  is  a  slight  improvement 
compared  to  the  theoretical  result  and  j|ACj||  decreases  slightly.  However,  as  Mach 
number  continues  to  increase,  a  pronounced  degredation  in  the  results  ensues,  find 
||AC/j|  increases  by  an  order  of  magnitude  (case  11).  Cases  12  through  15  examined 
the  trend  for  increasing  Reynolds  number  (beyond  3  x  106).  The  Reynolds  num¬ 
ber  was  increased  in  increments  of  3  x  106  up  to  15  x  106.  The  results  indicate  a 
continued,  albeit  gradual,  convergence  to  the  theoretical  result.  Finally,  in  cases  16 
through  19,  a  limited  grid  refinement  effort  was  undertaken.  A  slight  degredation 
in  accuracy  was  noted  for  the  finer  grids,  when  compared  with  the  coarse  grid,  as 
||AC/||  increased  slightly  (cases  17  and  19).  This  may  be  due,  in  part,  to  a  better 
resolution  of  the  wake  and  boundary  layer,  where  the  resolved  viscous  effects  may 
have  an  influence.  It  should  be  noted  that  for  grids  finer  than  60x25,  the  turbulence 
model  must  be  employed  to  prevent  secondary  oscillations  from  developing.  The 
secondary  oscillations  further  increase  the  ||AC/||  norm,  as  indicated  in  cases  16  and 
18.  The  development  of  the  secondary  oscillations,  due  principally  to  the  reduced 
levels  of  effective  dissipation  when  the  turbulence  model  is  not  applied,  is  shown  in 
Figure  37. 

Also  displayed  in  Figures  40-43,  at  r  =  1,2,  and  3  are  the  vorticity  contours  in 
the  wake.  The  solutions  were  obtained  with  the  incompressible  code  using  Grid  2  at 
Re  —  500  and  k  =  2,  starting  from  an  equilibrium  solution  at  r  =  0  and  a  =  0,  with 
a  peak  angular  amplitude  of  approximately  2  degrees.  Positive  vorticity  is  generated 
on  the  lower  surface,  while  negative  vorticity  predominates  near  the  upper  surface. 
The  effect  of  the  oscillations  on  the  wake  development  is  pronounced,  as  the  shed 
vorticity  indicates.  The  aforementioned  difficulty  concerning  rj-upwinding  is  evident 
in  the  wake.  That  is,  with  no  ^-upwinding  applied,  there  is  no  effective  damping  in 
the  rj  direction  except  for  that  provided  by  natural  dissipation.  A  close  examination 
of  the  contours  in  the  wake  reveals  negligible  oscillations  occuring  along  the  contours 
in  the  £  direction.  However,  increased  numerical  oscillations  are  noted  when  the 
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Case  No. 

gtgggEI 

Re 

Mach 

/« 

/. 

turb 

II  AC,  || 

1 

60x25 

0.20 

0.0200 

0.0400 

■a 

gawawBiriMi 

2 

60x25 

0.20 

0.0100 

0.0200 

off 

IPEPEMI 

3 

60x25 

flWl 

0.20 

0.0050 

0.0100 

off 

4 

60x25 

0.20 

0.0025 

0.0050 

off 

diverged 

5 

60x25 

0.20 

0.0200 

0.0400 

on 

6 

60x25 

w 

0.25 

0.0200 

0.0400 

off 

7 

60x25 

0.30 

0.0200 

0.0400 

off 

8 

60x25 

KQBfl 

0.35 

0.0200 

0.0400 

menmi 

9 

60x25 

Kwmum 

0.40 

0.0200 

0.0400 

off 

10 

60x25 

0.45 

0.0200 

0.0400 

off 

mnaEianiHai 

1  11 

60x25 

| 

0.50 

0.0200 

0.0400 

off 

BEgEBjHI 

12 

60x25 

uam 

0.20 

0.0200 

0.0400 

off 

13 

60x25 

0.20 

0.0200 

0.0400 

warn 

!  14 

60x25 

Bl^STia 

0.20 

0.0200 

0.0400 

15 

60x25 

0.20 

0.0200 

0.0400 

off 

16 

125x50 

0.20 

0.0200 

0.0400 

off 

17 

125x50 

wmm 

0.20 

0.0200 

0.0400 

on 

18 

208x108 

0.20 

0.0200 

0.0400 

off 

19 

208x108 

mam 

0.20 

0.0200 

0.0400 

on 

Table  3.  Tabulated  Parametric  Results;  Reduced  Frequency,  k  =  1 


contours  are  oriented  in  the  rf  direction.  The  oscillations  cannot  be  solely  attributed 
to  grid  courseness  in  the  wake  region  because,  for  the  particular  region  in  question, 
the  grid  spacing  in  the  £  direction  is  actually  coarser  than  that  in  the  r?  direction, 
particularly  near  the  branch-cut  (cf.  Figure  5).  A  comparison  with  the  O-Grid  case 
is  offered  in  Figure  43,  using  the  compressible  code  at  the  same  conditions.  In  general 
the  development  of  vorticity  near  the  trailing  edge  and  the  subsequent  convection 
into  the  wake  is  less  pronounced.  There  are  two  possible  reasons  for  this.  Primarily, 
the  differences  in  topology  are  considerable  at  the  trailing  edge.  The  O-grid  has  a 
rounded  trailing  edge  while  the  C-grid  ends  in  a  sharp  point.  The  sharp  trailing  edge 
in  the  case  of  the  C-grid  appears  to  act  as  a  vortex  generator  as  the  airfoil  oscillates. 
Another  significant  difference  is  the  application  of  the  boundary  conditions.  The 
branch  cut,  for  the  case  of  the  C-grid,  extends  from  the  trailing  edge  to  the  outer 
computational  domain.  In  contrast,  the  branch  cut  for  the  O-grid  extends  from  the 
leading  edge  to  the  computational  boundary. 
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Figure  27.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Compressible  Code,  a  (short-dashed  line),  Reduced  Fre¬ 
quency  k  =  0.05;  (a)  Re  =  100,  (b)  Re  =  1000,  (c)  Re  =  1.5  x  10s,  (d) 
Re  =  3.0  x  106 


Figure  28.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Compressible  Code,  a  (short-dashed  line),  Reduced  Fre¬ 
quency  k  =  0.2;  (a)  Re  =  100,  (b)  Re  =  1000,  (c)  Re  =  1.5  x  105,  (d) 
Re  =  3.0  x  10* 
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Figure  29.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Compressible  Code,  a  (short-dashed  line),  Reduced  Fre¬ 
quency  k  =  1;  (a)  Re  =  100,  (b)  Re  =  1000,  (c)  Re  =  1.5  x  105,  (d) 
Re  =  3.0  x  106 
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Figure  30.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Compressible  Code,  a  (short-dashed  line),  Reduced  Fre¬ 
quency  k  =  2;  (a)  Re  =  100,  (b)  Re  =  1000,  (c)  Re  =  1.5  x  105,  (d) 
Re  =  3.0  x  10® 
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Figure  32.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Compressible  Code,  a  (short-dashed  line),  Reduced  FYe- 
quency  k  =  6;  (a)  Re  =  100,  (b)  Re  =  1000,  (c)  Re  =  1.5  x  10s,  (d) 
Re  =  3.0  x  10* 
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Figure  33.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Incompressible  Code,  a  (short-dashed  line),  Reduced 
FYequency  k  =  0.05;  (a)  Re  =  100,  (b)  Re  =  500,  (c)  Re  =  1000,  (d) 
Re  =  2000 
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Figure  34.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Incompressible  Code,  a  (short- dashed  line),  Reduced 
Frequency  k  =  0.2;  (a)  Re  =  100,  (b)  Re  =  500,  (c)  Re  =  1000,  (d) 
Re  =  2000 
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Figure  35.  Comparison  of  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Incompressible  Code,  a  (short-dashed  line),  Reduced 
Frequency  k  =  2;  (a)  Re  =  100,  (b)  Re  =  500,  (c)  Re  =  1000,  (d) 
Re  =  2000 


87 


5.0 

Figure  36. 
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Comparison  ;  Calculated  (solid  line)  and  Theoretical  (long-dashed  line) 
Lift  Coefficient,  Incompressible  Code,  a  (short-dashed  line),  Reduced 
Frequency  k  =  4;  (a)  Re  =  100,  (b)  Re  =  500,  (c)  Re  =  1000,  (d) 
Re  =  2000 
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Figure  37.  O-Grid  209x108:  k  =  1,  Re  =  3  x  10®,  (a)  Turbulence  Model  Off,  (b) 
Turbulence  Model  On 
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Figure  38.  Phase  and  Amplitude  Errors  as  a  Function  of  Reduced  Frequncy:  Com¬ 
pressible  Code;  Phase  Error  (O)  Amplitude  Error  (A) 
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Figure  43.  Vorticity  Contours:  Compressible  Code,  O-grid  120x50,  Re  =  500, 
k  =  2,  r  =  3.0 

5.4  Conclusions 

The  influence  of  the  turbulence  model  is  demonstrated  in  the  numerical  results 
for  the  incompressible  code.  Implementation  of  the  modified  Baldwin- Lomax  turbu¬ 
lence  model  allowed  the  numerical  calculation  of  flutter-onset  speeds  corresponding 
to  theoretical  aud  experimental  values  for  a  NACA  0015  airfoil.  The  numerical  on¬ 
set  of  flutter  is  a  function  of  u)k  and  u>a  independently,  in  contrast  to  the  theoretical 
results,  which  predict  onset  as  a  function  of  the  ratio  of  these  parameters.  The  rea¬ 
son  for  this  effect  is  the  subcritical  nature  of  the  bifurcation.  There  is  experimental 
evidence  of  this  as  well,  where  Yang  and  Zhao  [11]  noted  a  plunge-dominated  mode 
that  could  only  be  excited  by  a  physical  displacement  of  the  airfoil,  indicating  a 
differential  response  to  perturbation  level. 

The  greatest  challenge  to  the  accurate  calculation  of  flutter  onset  for  a  given 
structure  is  the  development  of  turbulence  models  applicable  in  this  regime.  The 
Baldwin-Lomax  model  has  been  criticized  in  the  unsteady  and  complex  flow  regimes 
[14, 47],  and  in  the  abilit”  to  accurately  predict  self-induced  oscillations  [16].  It  is  un¬ 
clear  at  present  whether  alternatives,  such  as  the  Jfc  —  c  model  [47],  would  be  superior. 
The  present  research  indicates  that  the  Baldwin-Lomax  model  may  serve  adequately 
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in  predicting  flutter  onset  associated  with  bifurcations  from  an  equilibrium  state.  In 
such  a  regime,  the  flow  is  not  dynamic  and  the  application  of  the  turbulence  model 
may  be  more  appropriate.  This  conclusion  is  substantiated  numerically  for  the  case 
of  forced,  small  amplitude  airfoil  oscillations  where  the  inclusion  of  the  turbulence 
model  produces  no  discemable  impact  upon  the  computed  lift  time  history.  How¬ 
ever,  since  the  scope  of  this  demonstration  is  limited,  further  research  is  warranted 
in  this  area. 

A  basis  is  established  to  validate  the  aerodynamic  response  for  the  compressible 
and  incompressible  codes  and  their  application  to  the  numerical  calculation  of  flut¬ 
ter.  The  strong  correspondence  between  the  lift  coefficient  predicted  by  Theodorsen’s 
function  and  the  numerically  calculated  lift  coefficient  across  a  range  of  reduced  fre¬ 
quencies  indicates  the  applicability  of  these  codes  in  this  capacity.  The  incompress¬ 
ible  code  is  superior  to  the  compressible  code  in  the  low  Reynolds  number  regime, 
while  the  compressible  code  is  superior  across  the  range  of  frequencies  and  at  higher 
Reynolds  numbers.  The  incompressible  code  fails  to  accurately  predict  the  unsteady 
lift  for  k  >  4.  This  failure  is  attributed  to  the  lack  of  artificial  dissipation  in  the 
T)  direction.  The  inability  to  apply  simultaneous  if-  and  £-upwinding  for  the  non- 
conservative  form  of  the  equations  is  therefore  a  detriment.  Possible  solutions  to  this 
problem  are  rewriting  the  equations  in  conservative  form,  the  application  of  artificial 
dissipation  in  conjuction  with  central-differencing,  or  both.  The  compressible  code 
also  displays  increased  error  associated  with  increased  frequencies,  k  >  6,  but  they 
are  less  pronounced. 
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VI.  Summary 


This  section  is  presented  to  provide  a  synthesis  of  the  results  and  conclusions 
provided  in  Parts  I  and  II. 

In  the  laminar  flow  regime  explored  in  Part  I,  the  influence  of  the  structural 
model  of  the  airfoil  on  the  location  of  the  Hopf  bifurcation  structure  has  been  demon¬ 
strated.  In  general,  the  structural  model  allows  the  development  of  sustained,  limit- 
cycle  oscillations  in  the  flowfield  to  develop  at  a  lower  Reynolds  number  than  would 
be  the  case  for  a  fixed  airfoil.  This  effect  is  most  evident  when  the  structural  model 
incorporates  multiple  degrees  of  freedom.  Coupling  between  the  axes  can  occur 
through  either  direct  structural  coupling,  or  indirectly  via  the  aerodynamics.  Single¬ 
axis  experiments  demonstrate  that  the  vertical  axis  has  a  minimal  impact  on  the 
value  of  Reynolds  number  at  which  the  flow  becomes  oscillatory  in  the  absense  of 
a  resonance  effect,  while  the  pitch  axis  has  a  more  pronounced  effect.  The  effect  of 
structural  damping  is  most  pronounced  in  the  bistable  region,  beyond  this  region,  the 
structural  damping  impacts  the  airfoil  motion  more  strongly  than  the  aerodymanic 
flowfield,  which  will  remain  unsteady  regardless  of  the  structural  damping  level. 

The  structure  of  the  bifurcation  is  subcritical,  as  the  existence  of  perturbations 
in  the  flowfield  can  initiate  sustained  oscillations  at  a  Reynolds  number  lower  than 
would  occur  when  starting  from  an  initial  equilibrium  state.  The  influence  of  the 
structural  model  may  be  more  pronounced  on  the  branch  of  periodic  solutions  in  the 
bistable  region  than  on  the  Hopf  point  itself.  In  addition,  the  structural  model  is  seen 
to  provoke  additional  frequency  content  in  the  time  history  of  C/,  for  example.  The 
addition  of  structural  elements  may  induce  secondary  bifurcations  which  perturb  the 
solution  trajectory  from  a  single,  closed  orbit  in  phase  space  to  a  more  complicated 
orbit.  This  is  true  in  the  long-term  time  behavior  and  in  the  immediate  evolution 
away  from  an  equilibrium  solution. 
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The  influence  of  the  turbulence  model  is  demonstrated  in  the  numerical  results 
for  the  incompressible  code  in  Part  II.  Implementation  of  the  modified  Baldwin- 
Lomax  turbulence  model  allowed  the  numerical  calculation  of  flutter-onset  speeds 
corresponding  to  theoretical  and  experimental  values  for  a  NACA  0015  airfoil.  The 
numerical  onset  of  flutter  is  a  function  of  and  u>a  independently,  in  contrast 
to  the  theoretical  results,  which  predict  onset  as  a  function  of  the  ratio  of  these 
paratmeters.  The  reason  for  this  effect  is  the  subcritical  nature  of  the  bifurcation. 
There  is  experimental  evidence  of  this  as  well,  where  Yamg  and  Zhao  [11]  noted  a 
plunge-dominated  mode  that  could  only  be  excited  by  a  physical  displacement  of  the 
aurfoil,  indicating  a  differential  response  to  perturbation  level. 

The  greatest  challenge  to  the  accurate  calculation  of  flutter  onset  for  a  given 
structure  is  the  development  of  turbulence  models  applicable  in  this  regime.  The 
Baldwin-Lomax  model  has  been  criticized  in  the  unsteady  and  complex  flow  regimes 
[14,  47],  amd  in  the  ability  to  accurately  predict  self- induced  oscillations  [16].  It  is 
unclear  at  present  whether  alternative  turbulence  models,  such  as  the  k—e  model  [47], 
would  be  superior.  The  present  research  indicates  that  the  Baldwin-Lomax  model 
may  serve  adequately  in  predicting  flutter  onset  associated  with  bifurcations  from 
an  equilibrium  state.  In  such  a  regime,  the  flow  is  not  dynamic  amd  the  application 
of  the  turbulence  model  may  be  more  appropriate.  This  conclusion  is  substantiated 
numerically  for  the  case  of  forced,  small  amplitude  airfoil  oscillations  where  the 
inclusion  of  the  turbulence  model  produces  no  discemable  impact  upon  the  computed 
lift  time  history.  However,  since  the  scope  of  this  demonstration  is  limited,  further 
research  is  warramted  in  this  area,  particularly  in  the  area  of  grid  refinement,  as  it 
was  noted  on  the  finer  grids  that  application  of  turbulence  inhibits  the  bifurcation 
to  secondary  oscillations.  The  fundamental  mode  of  oscillation  remains  unchamged 
for  both  coarse  amd  fine  grids. 

A  basis  is  established  to  validate  the  aerodynamic  response  for  the  compressible 
and  incompressible  codes  and  their  application  to  the  numerical  calculation  of  flut- 
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ter.  The  strong  correspondence  between  the  lift  coefficient  predicted  by  Theodorsen’s 
function  and  the  numerically  calculated  lift  coefficient  across  a  range  of  reduced  fre¬ 
quencies  indicates  the  applicability  of  these  codes  in  this  capacity.  The  incompress¬ 
ible  code  is  superior  to  the  compressible  code  in  the  low  Reynolds  number  regime, 
while  the  compressible  code  is  superior  across  the  range  of  frequencies  and  at  higher 
Reynolds  numbers.  The  incompressible  code  fails  to  accurately  predict  the  unsteady 
lift  for  k  >  4.  This  failure  is  attributed  to  the  lack  of  artificial  dissipation  in  the 
rf  direction.  The  inability  to  apply  simultaneous  r?-  and  £-upwinding  for  the  non¬ 
conservative  form  of  the  equations  is  therefore  a  detriment.  Possible  solutions  to  this 
problem  are  rewriting  the  equations  in  conservative  form,  the  application  of  artificial 
dissipation  in  conjuction  with  central-differencing,  or  both.  The  compressible  code 
also  displays  increased  error  (when  compared  with  the  theoretical  result)  associated 
with  increased  frequencies,  k  >  6,  but  they  are  less  pronounced. 

Grid  sensitivities  to  critical  parameters  such  as  the  Strouhal  number  and  crit¬ 
ical  Reynolds  number  are  considerable.  Relatively  minor  changes  in  grid  spacing 
or  domain  size  may  have  a  sizeable  effect  on  numerical  results.  At  a  minimum,  a 
limited  grid  refinement  study  is  therefore  required  to  establish  consistent  results. 
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Appendix  A.  Code  Validation 


The  validation  of  the  incompressible  code  was  divided  into  two  parts.  In  the 
first  part,  as  reported  by  Lutton  and  Beran  [5],  the  accuracy  of  the  unsteady  nu¬ 
merical  results  was  examined  for  the  case  of  a  circular  cylinder.  The  cylinder  was 
chosen  because  of  the  availability  of  experimental  and  numerical  results.  In  the  sec¬ 
ond  part  of  the  validation,  the  application  of  the  vertical  axis  was  examined  within 
the  context  of  the  aerodynamic  portion  of  the  code.  The  vertical  axis  was  examined 
singularly  because  the  equation  for  the  pitch  axis  is  of  identical  form. 

A.l  Integration  of  the  Navier-Stokes  Equations 

The  algorithm  and  the  implementation  of  the  algorithm  are  validated  by  apply¬ 
ing  the  numerical  procedure  to  the  simulation  of  flow  about  a  fixed  circular  cylinder. 
Details  specific  to  the  application  are  discussed  first,  followed  by  a  comparison  of 
computed  results  with  data  reported  in  previous  experimental  and  numerical  inves¬ 
tigations. 

The  grid  structure  used  for  the  analysis  of  an  airfoil  flow  was  not  appropri¬ 
ate  for  a  similar  analysis  involving  a  circular  cylinder,  since  a  sharp  trailing  edge  is 
absent  from  the  latter  geometry.  For  the  circular  cylinder,  the  numerical  procedure 
was  modified  to  treat  an  O-grid,  like  the  one  pictured  in  Figure  44a.  Aside  from 
trivial  changes  in  the  numerical  algorithm  to  account  for  the  altered  grid,  substantial 
changes  in  the  application  of  numerical  dissipation  and  boundary  condition  enforce¬ 
ment  were  required.  Also,  an  alternative  grid  generation  technique  was  employed. 
In  the  development  of  the  governing  equations  for  the  new  flowfield,  the  radius  of 
the  cylinder,  a,  was  chosen  as  a  length  scale.  With  the  velocity  field  scaled  by  U , 
a  natural  Reynolds  number  was  obtained,  Re„  =  Ua/v ,  which  is  proportional  to  a 


97 


Figure  44.  O-grid  Structure:  (a)  Node  Distribution;  (b)  Schematic  of  Boundaries 
Reynolds  number  based  on  cylinder  diameter,  d: 

Ud 

Re*  =  —  =  2Re0. 
v 

For  the  O-grid,  node  points  were  distributed  along  rays  emanating  from  the 
cylinder  center.  The  rays  were  clustered  about  the  grid  cut,  AB,  aft  of  the  cylinder 
(see  Figure  44b).  In  an  equivalent  manner  on  each  ray,  the  radial  node  spacing  was 
increased  through  a  geometric  progression,  starting  with  a  minimum  value,  As„, 
at  the  cylinder  surface.  The  radial  position  of  the  outer  boundary,  R,  the  number 
of  nodes  in  the  radial  direction,  J,  and  As„,  serve  to  define  uniquely  the  node 
arrangement.  The  angular  distribution  of  rays  was  also  calculated  using  a  geometric 
progression,  where  the  smallest  angular  increment,  2/?x/(  J- 1),  was  specified  to  occur 
adjacent  to  the  cut,  while  the  largest  increment  was  specified  to  occur  adjacent  to 


the  forward  stagnation  point  (nodes  were  placed  symmetrically  about  the  centerline). 
Here,  I  is  the  number  of  nodes  in  the  azimuthal  direction  and  $  is  &  free  parameter 
that  was  chosen  to  be  0.5. 

Results  were  obtained  using  three  grids,  01  (coarse  grid),  02  (fine  grid)  and 
03  (large  and  fine  grid).  Characteristics  of  these  grids  are  summarized  in  Table  4. 

The  O-grid  structure  differs  most  from  the  C-grid  structure  in  the  wake  region. 
There,  lines  of  constant  17  are  basically  aligned  with  the  freestream  direction  for 
the  airfoil  and  normal  to  the  freestream  for  the  circular  cylinder.  Owing  to  the 
importance  of  numerical  dissipation  in  the  computation  of  smooth  wake  flows,  it 
was  necessary  to  switch  the  application  of  upwinding  from  £  derivatives  (airfoil)  to 
77  derivatives  (circular  cylinder).  The  change  in  the  application  of  upwinding  was 
straightforward  to  implement,  since  the  convection  terms  are  handled  explicitly. 

Outflow  conditions  in  the  airfoil  analysis  were  naturally  enforced  along  the 
straight  boundary  of  the  C-grid,  downstream  of  the  airfoil.  For  the  O-grid,  outflow 
conditions  were  specified  over  a  predefined  arc,  CD  (see  Figure  44b).  The  number 
of  nodes  comprising  the  arc  was  2/7,  where  for  all  calculations  7  =  0.15  was  chosen. 
The  vorticity  condition  enforced  on  arc  CD  was  =  0,  which  was  evaluated  using  a 
2-point  approximation  to  u)v.  Sa  and  Chang  [50]  applied  the  same  outflow  condition, 
but  over  an  arc  defined  by  the  local  properties  of  rp. 

After  the  integration  algorithm  was  modified  to  incorporate  an  O-grid  struc¬ 
ture,  flows  over  a  fixed,  circular  cylinder  were  numerically  simulated  for  Reynolds 
numbers  between  10  and  80  in  increments  of  10.  Drag  coefficient  and  Strouhal  num¬ 
ber,  St  =  f,d/U,  were  used  to  compare  computed  results  to  data  available  in  the 


Grid 

R 

mm 

mm 

bs 

01 

101 

51 

51 

0.01 

02 

101 

101 

101 

0.01 

03 

201 

101 

201 

0.01 

Table  4.  Characterstics  of  Grids  Ol,  02  and  03 
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Figure  45.  Strouhal  Number  u  a  Function  of  Reynolds  no.  Experimental:  Roshko 
[28],  Tritton  [29]  and  Berger  [30];  Computational:  Sa  and  Chang  [51], 
and  Present  (Grids  02  and  03). 

literature.  Time-periodic  solutions  were  obtained  at  Reynolds  numbers  equal  to  or 
greater  than  50.  No  attempt  was  made  to  compute  the  critical  Reynolds  number  (be¬ 
tween  40  and  50)  at  which  the  flow  first  becomes  unsteady.  A  time  step  of  0.04  was 
used  for  Reynolds  numbers  greater  than  50.  At  smaller  Reynolds  numbers,  for  which 
the  computed  flowfields  were  steady,  larger  time  steps  were  taken.  Time-dependent 
results  were  insensitive  to  reductions  of  time  step  below  0.04.  Solutions  were  also 
found  to  be  insensitive  to  changes  in  the  geometric  parameter  As„.  Sensitivity  to 
the  parameters  and  7  was  not  evaluated. 

Computed  Strouhal  number  was  compared  to  the  results  of  experimental  and 
numerical  investigations  in  Figure  45.  It  should  be  noted  that  the  solutions  obtained 
with  grid  01  for  Re,*  >  50  were  unsteady,  but  not  time-periodic,  and  behaved  in  a 
spurious  manner.  It  is  believed  that  this  departure  from  periodicity,  not  observed 
with  grids  02  and  03  for  the  range  of  Reynolds  number  reported,  is  directly  at- 
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tributable  to  grid  coarseness  and  the  non-conservative  formulation  of  the  convective 
terms.  Thus,  aperiodic  solutions  obtained  with  grid  01  were  not  included  in  the 
results  presented  below.  The  experimental  studies  of  Roshko  [28],  Tritton  [29]  and 
Berger  [30]  and  the  numerical  investigation  of  Sa  and  Chang  [51]  served  as  a  basis 
for  comparison.  The  computed  results  were  found  to  be  in  excellent  agreement  with 
Roshko’s  data,  except  for  grid  02  near  Rej  =  50.  Apparently,  at  Reynolds  numbers 
just  exceeding  the  critical  Reynolds  number,  the  effects  of  domain  size  are  signifi¬ 
cant,  and  a  domain  size  of  R  =  101  is  insufficient.  The  trend  of  decreasing  St  with 
increasing  R  was  also  observed  by  Sa  and  Chang  [50]  when  the  freestream  boundary 
condition  on  perturbation  streamfunction,  tp  =  0,  was  enforced.  It  is  unclear  to  the 
authors  whether  the  disparity  between  the  present  data  and  that  computed  by  Sa 
and  Chang  [51]  is  a  result  of  differences  in  the  treatment  of  the  convective  terms  or 
in  differences  of  grid  structure.  Sa  and  Chang  employ  a  conservative  formulation 
of  the  convective  terms,  with  4th-order  approximations  to  velocity  components,  and 
report  results  for  a  51  x  50  O-grid. 

The  drag  coefficients  associated  with  computed  flowfields  about  a  circular  cylin¬ 
der  are  plotted  in  Figure  46  versus  Reynolds  number.  The  experimental  study  of 
Tritton  [29]  and  the  numerical  investigations  of  Borthwick  [52]  and  Sa  and  Chang 
[51]  were  used  as  a  basis  for  comparison.  When  Rej  <  40  (steady  flow),  the  results 
of  the  present  study  show  little  sensitivity  to  grid  parameters.  Computed  values  of 
Cd  are  5%  to  10%  below  Tritton’s  best-fit  data.  At  a  Reynolds  number  of  50  there 
is  a  noticeable  difference  between  data  obtained  with  grids  02  and  03.  The  larger 
domain,  03,  provides  a  close  match  to  the  value  predicted  by  Sa  and  Chang,  perhaps 
owing  to  their  enforcement  of  a  highly  accurate  far-field  condition  on  streamfunction. 
As  Reynolds  number  was  increased  beyond  50,  computed  drag  values  obtained  with 
grids  02  and  03  were  quite  close,  but  tended  to  somewhat  exceed  that  indicated 
by  experiment.  Numerical  integration  at  Reynolds  numbers  beyond  about  150  is 
complicated  by  the  development  of  turbulence  in  the  cylinder  wake  (Roshko  [28]). 
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Figure  46.  Drag  Coefficient  as  a  Function  of  Reynolds  Number.  Experimental: 

Tritton  [29];  Computational:  Sa  and  Chang  [51],  Borthwick  [52],  and 
Present  (Grids  01,  02  and  03). 
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Figure  47.  Comparison  of  Exact  (solid  line)  and  Runge-Kutta  (©)  Solutions:  Ver¬ 
tical  Axis 


A. 2  Application  of  Structural  Model 

As  mentioned  earlier  (cf.  Section  2.5),  the  aerodynamic  coefficients  are  calcu¬ 
lated  immediately  after  the  calculation  of  and  w,  and  are  treated  as  constants  over 
the  interval  t  to  t  +  At.  This  introduces  a  phase  lag  of  up  to  At  in  the  application 
of  the  loads.  However,  at  the  small  timesteps  required  for  the  stability  of  the  overall 
scheme,  this  does  not  introduce  significant  error.  To  substantiate  this  statement,  the 
vertical  axis  (cf.  Eq.  (18))  was  examined  to  validate  the  Runge-Kutta  integration 
procedure.  A  prescribed  forcing  function  was  applied  to  the  right-hand  side  of  Eq. 
(18): 

Qh  =  A0  sin(u;/t  +  6).  (143) 

This  equation  has  an  exact  solution  (for  6  =  0)  given  by  [53] 

M‘)=;d^ain(“''+9)'  (144) 

where  a  =  2 £r,  6  =  1  —  ra,  r  =  vj/uh,  and  6  =  arctan(u/6).  Values  for  the  constants 
were  m0  =  1.0,  Sa  =  0.0,  Dh  =  20.0,  and  Kh  =  25.0,  for  which  £  =  2.0.  The  solu- 
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tion  was  then  integrated  numerically  by  applying  the  Runge-Kutta  method  using  the 
same  forcing  function  except  that  in  this  case  8  was  set  to  0.01.  This  value  of  8  repre¬ 
sented  the  largest  timestep  employed  throughout  the  numerical  experiments.  Most 
timesteps  were  considerably  smaller,  particularly  when  finer  grids  were  employed. 
The  results  are  displayed  in  Figure  47.  No  pronounced  difference  is  evident  in  the 
two  solutions,  and  it  is  judged  that  the  application  of  the  lagged  forcing  function 
with  Runge-Kutta  integration  has  a  minimal  impact  on  airfoil  motion. 
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Appendix  B.  Comparison  of  Incompressible  Codes 

B.l  Validity  of  LU  Decomposition  Approach 

The  method  of  solution  for  the  linear  system  of  equations  presented  in  Chap¬ 
ter  2  (Eq.  (59)) 

AAntl>i  =  N4  +  Gi(L2N2  -  Nr)  - 

wb~re 

A  =  -(GiLi  -f-  GtL3), 

was  chosen  to  be  LU  decomposition  because  the  matrix  A  is  time-invariant  (for 
a  constant  time  step  and  Reynolds  number),  and  can  therefore  decomposed  into 
the  product  of  lower  and  upper  triangular  matrices  once  at  the  start  of  the  time- 
integration  procedure. 

The  operators  Li  and  L2  represent  the  discrete  form  of  the  Laplace  operator 

A 

acting  on  Arp  to  satisfy  the  Poisson  equation.  Further,  the  operator  G  implicitly 
applies  the  Laplace  operator  to  Anu>.  However,  the  Poisson  equation  (Eq.  (28)), 
shown  here  again  for  convenience, 

J2V2rp  =  fatpx  +  <pxrpm  -  2 fcrpiv  +  <p4rp£  +  <f>srp„  =  -j2v, 

contains  the  terms  <f>\  through  fa: 

<Pi  =  x\  +  y£, 

<h  =  *;+ij, 

<h  =  x£x„  +  y£yv, 

<t>*  =  */_1  [My™**  -  xmy*)  +  Mvuxn  ~  *«<y.>)  +  2  6»(*€„y„  - 

fa  =  -3~l [^i(y*mxt  ~  xwVi)  +  Myttxt  -  xuy()  +  2<h(xtr,yt  -  y«»»x€ )] - 
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It  is  then  required  to  demonstrate  the  validity  of  the  described  approach  when  the 
time-dependent  metric  transforation  introduced  in  Chapter  4  (Eq.  (69))  is  applied. 
In  this  case  the  metrics  are  functions  of  time,  and  the  <j>i  are  functions  of  the  metrics 
and  their  derivatives.  To  this  end,  consider  first  the  Jacobian  J 

J  =  X{yv  -  x„yt  #  0, 

since  it  is  intuitively  clear  that  for  the  case  of  a  rigid  grid,  the  Jacobian  must  remain 
time-invariant  under  the  action  of  arbitrary  grid  motions.  Imposing  a  prescribed 
coordinate  transformation  from  system  (x,  y)  to  system  (x,y): 

x  =  x  cos  0  —  y  sin  0  -f  hx , 
y  =  x  sin  0  +  y  cos  0  +  hyy 

where  0  represents  the  angle  of  rotation  and  hx,  hy  the  respective  translational  dis¬ 
placements.  In  this  case,  the  metrics  are  related  by 

cos  0  —  yt  sin  0,  xv  =  xv  cos  0  —  y„  sin  0, 
y*  =  !(_  sin  0  +  y*  cos  0 ,  y„  =  xv  sin  0  +  yv  sin  0 , 

and  further 

=  %  cos  0  -  y#  sin  0,  x„v 

y«  =  xk  sin  0  -f  y{f  cos  0 ,  y„„ 

X(v  =  z«*7  cos  0  -  yiv  sin  0 ,  y*„ 

Substituting  these  into  the  definition  of  J 

J  =  (x*co8  0  —  y^  sin  0)(xnsin0  +  sir.  f'< 


=  xvv  cos  0  -  yvv  sin  0 , 
=  x^  sin  0  -I-  ym  sin  0 , 
=  x^sinfl  +  y^r,  sin  0, 
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—  (fa  cos0  —  y,sind)(if,sind  +  y,,sin0) 

=  if y, (cos2  0  +  sin2  0)  -  i,yf  (cos2  0  +  sin2  8) 

+  if  i,  cos  0  sin  0  —  jjfy„cos0sin0 

—  if  i,  cos  0  sin  0  +  yf  y„  cos  0  sin  0 

=  fayv-faytt 

and  so  the  Jacobian  is  time-invariant,  as  expected.  Now  the  process  is  repeated  for 
the  fa ,  starting  with  fa: 

fa  =  (if  cos  0  —  y f  sin  0)2  +  (if  sin  0  +  yf  cos  0)2 
=  i|  (cos2  0  +  sin2  0)  +  y|  (cos2  0  +  sin2  0) 

+  2if  yf  cos  0  sin  0  —  2if  yf  cos  0  sin  0 

=  *f  +  y?, 


and  continuing  with  fa 

fa  =  (i„cos0  -ynsin0)2  +  (i„  sin  0  +  y,,  cos  0)2 
=  i2(cos2  0  +  sin2  0)  +  y2(cos2  0  +  sin2  0) 

+  2xvyv  cos  0  sin  0  -  2iTJyl,  cos  0  sin  0 

=  *?  +  yJ, 

and  fa 

fa  =  (if  cos0  -  yf  sin0)(i,,cos0  -  y,sin0) 

+  (if  sin  0  +  yf  cos  8)(fa  sin  0  -f  yn  cos  0) 

=  ifi„(cos2  0  +  sin2  0)  +  yfy,(cos2  0  +  sin2  0) 
+  if  y„  cos  0  sin  0  —  ify„cos0sin0 
-I-  i„yf  cos  0  sin  0  —  i,,yf  cos  0  sin  0 
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=  fafa  +  tkiiv 


With  J,fa,fa  and  fa  established  as  time-  invariant,  all  that  remains  to  demonstrate 
the  same  for  fa  and  fa  is  to  examine  the  terms  containing  second  derivatives  in  their 
definition,  that  is, 


xrmVv) 


(xw  sin  0  +  ym  cos  0)(xv  cos  0  -  y„  sin  0) 
(xw  cos  0  -  y„„  sin  0)(x„  sin  0  +  yv  cos  0) 
xvnXr,  cos  0  sin  0  +  y^x,  cos2  0 
xnvyv  sin2  0  -  yrnijr,  cos  0  sin  0 


+  y,„x„sin20  -  x,„x„cos0sm0 

-  x~r,yv  cos2  0  +  ymyn  cos  0  sin  0 

—  (y^fa  xwVv)y 


and  by  permutation  of  subscripts,  it  is  established  in  similar  fashion  that 


-  *«y„) 

(XivVri  ~  y^vXv) 
iVrmxi  ~  xv'iVl ) 
(y«X{  -  x^yc) 

{hm  -  y«„*<) 


(faifa  -  *«y»), 

(far,  y„  -  favfah 

(yr>r>X(  ~  Xrr?y()i 

(mfa  -  faifa)r 
(xinfa  ~  tiivfa)- 


It  is  therefore  established  that  the  application  of  LU  decomposition  as  described  is 
valid. 


B.2  Stability  Analysis  by  Model  Equation 

To  investigate  the  difference  in  the  stability  characteristics  of  the  two  versions 
of  the  incompressible  code,  a  comparative  stability  analysis  is  undertaken.  A  com- 
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parative  stability  analysis  of  this  type  is  not  available  in  the  literature.  The  stability 
analysis  is  accomplished  using  a  Fourier  analysis  (also  refered  to  as  the  von  Neu¬ 
mann  analysis)  [54]  applied  to  the  linearized,  viscous  Burger’s  equation.  Burger’s 
equation  is  selected  because  of  the  similarity  to  the  vorticity-transport  equation  (cf. 
Eq.  (29)). 

To  facilitate  the  examination  of  the  two  methods,  a  similar  treatment  is  applied 
to  Burger’s  equation.  In  the  first  case,  a  source  term  is  applied  which  represents 
the  apparent  body  forces  which  are  similarly  applied  to  the  momentum  equation 
(see  Chapter  2).  In  the  second  case,  a  time-dependent  coordinate  transformation 
is  applied  to  Burger’s  equation,  as  in  the  modified  form  of  the  incompressible  code 
(see  Chapter  4).  The  Fourier  analysis  is  applied  after  discretization  utilizing,  in  both 
cases,  a  forward-time  (first-order)  and  central-space  (second-order)  approach. 

Proceeding  with  the  Case  1,  Burger’s  equation  with  source  term 

S  =  yil- f  xft? 


with  y  representing  a  free  parameter,  is  written  as 


du  du  d2u 

s?  +  C5J  = ''ai5  +  • 


(145) 


The  discretized  equation  is  written  as 


u"+1  —  u? 

J  +c— 
At 


u”+i  -  u"-i 


2Ax 


Ax2 


or  upon  rearranging 


un+]  =  +  0(„?+1  _  2u”  +  +  S(A(). 


The  exact  solution  and  the  error  must  assume  the  same  form  [54],  hence  the  error  is 
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written  as  u” — ti^ac*)".  Assuming  an  error  distribution  of  the  form  em(x,  t)  =  e°*e,kmX 
where  m  is  a  summation  index  in  the  Fourier  series  representing  the  wavenumber, 
i.e., 

e(x,f)  =  Y2z™e'kmT' 

m 

with  zm  —*■  eat ,  we  have 


ea(t+At)eikx 


(eateik(*+Az) 
Ax 2  v 


j  cA<  /  ate«'t(x+A*)  _  eote«*(ar-Ax) 

2  Ax  ' 

e.tc.*(«-Ax))  + 


Dividing  the  above  equation  by  eatetkx 


and  using  the  identities 


cos/9  =  §(e*  +  e  ,/J)> 
t  sin  /?  =  5(e*/}  —  e~'p). 


where  j8  =  kTOx,  the  resulting  equation  can  be  written  as 


e 


oA< 


(‘  ~  20)  ~  '1^  “  ^  +  20  +  S(A‘> 
1  +  2r2(cos  y9  —  1)  —  *>i  sin  P  -f-  S(A<), 


with  rj  =  cAt/Ax  and  r2  =  fiAt/Ax 2. 

The  modulus  of  the  amplification  factor,  eaAt,  must  remain  less  than  unity, 
resulting  in  the  stability  restriction 


A\  =  |1  +  2r2(cos  P  —  1)  —  in  sin  P  -f  S(A<)|  <  1.  (146) 


A  similar  analysis  is  undertaken  for  Case  2.  The  time-dependent  coordinate 
transformation 


{  =  {(*,*)>  r  =  t, 
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is  applied  to  Eq.  (145)  in  lieu  of  the  source  term  applied  in  Case  1.  The  transformed 
equation  is  written  as 


8U 


du\ 

dt)' 


with  transformation  metrics  given  by 


and  the  Jacobian  defined  by 

J  = 

xi 

To  retain  the  same  form  for  the  viscous  terms,  choose  a  transformation  such  that 
Zx  =  1  and  Zxx  =  0>  resulting  in 

du  .du  _  d2u 
dT+Cd(~  ^ ' 

where  c  =  c  +  Zt  is  the  modified  convective  coefficient.  Note  that  the  redefinition  of 
the  convective  term  does  not  directly  correspond  to  the  inclusion  of  a  source  term. 

Application  of  the  Fourier  analysis  for  this  case  results  in  the  amplification 
factor  for  Case  2,  A2, 


A%  =  eoA<  =  1  +  2r2(cos  /?  —  1)  —  t(ri  +  — )  sin/?. 

Ax 


(148) 


Corresponding  to  the  case  of  pure  rotation,  set  xT  =  fly,  for  which  Zt  =  —  fty/J. 
In  comparing  the  amplification  factors  for  the  two  cases  (Eqs.  (146)  and  (148)),  it  is 
noted  that  both  retain  a  dependence  on  the  parameter  ft,  however,  the  source  term 
in  the  first  case  is  proportional  to  ft3.  In  addition,  the  explicit  dependence  on  ft 
vanishes  in  the  second  case.  It  is  also  noted  that  the  source  term  5  contributes  to 
the  red  part  of  At,  while  the  additional  terms  in  the  second  case  contribute  to  the 
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imaginary  part  of  A 2. 

The  parameters  for  the  two  cases  were  selected  to  represent  a  typical  appli¬ 
cation,  with  T\  =  0.5,  r2  =  0.375,  A t  —  .01  and  SI  =  ft  =  1.0.  In  Case  1  y  =  2.0 
corresponding  to  the  application  of  the  oscillatory  boundary  conditions  on  the  outer 
computational  boundary,  while  in  the  Case  2,  y  =  0.2  corresponding  with  the  ap¬ 
plication  of  the  oscillatory  boundary  conditions  on  the  airfoil  surface.  The  real  and 
imaginary  components  of  the  amplification  factors  Ai  =  eoAt  are  plotted  as  the  phase 
angle  /?  is  rotated  through  2ir.  The  results  are  shown  in  Figure  47.  The  dashed  line 
represents  the  unit  circle,  and  thus  represents  a  stability  boundary.  The  solid  line 
in  Figure  47a  indicates  the  modulus  of  A\.  The  results  show  that  a  growth  of  the 
amplification  factor  occurs  for  Case  1,  as  the  source  term  contributes  to  the  increase 
in  the  real  component.  In  the  Case  2,  however,  shown  in  Figure  47b  the  rotational 
rates  contribute  to  an  increase  in  the  imaginary  component,  where  there  is  a  greater 
stability  margin. 
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Appendix  C.  Direct  Calculation  of  Hopf  Points 


The  conditions  for  Hopf  bifurcation  are  given  by  Seydel  [6].  Let  F(y,  0)  repre¬ 
sent  a  non-linear  mapping  from  R?  x  R  into  i?",  with  y  €  R?  aad  ft  £  R.  A  solution 
of  the  system  of  differential  equations 

yt  +  F(y,P)  =  o 

is  steady-state  at  a  point  (y°,  f3°)  if 

F(y°,0°)  =  0. 


If  further  the  Jacobian  matrix 


A  = 


dF 

dy 


has  a  simple  pair  of  purely  imaginary  eigenvalues  A  =  ±i0,  and  the  so-called 
“transversality  condition”  is  satisfied 


±(Real{\m)  #0, 
dP  g. 

then  there  emerges  a  branch  of  limit-cycle  solutions.  The  period  of  the  limit  cycle  is 
given  by  2 ir/0,  and  the  initial  amplitude  is  zero. 

The  following  sections  detail  the  calculation  of  Hopf  bifurcation  points  for  the 
fixed  airfoil  [9]  and  the  modifications  made  for  the  case  of  the  moving  airfoil. 


C.l  Hopf  Bifurcation  Algorithm:  Fixed  Airfoil 

The  following  analysis  by  Beran  and  Lutton  [9]  permits  the  application  of  the 
method  of  Griewank  and  Reddien  [55]  to  the  direct  calculation  of  the  Hopf  point. 
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The  equations  of  motion  are  recast  by  using  an  operator  notation  (per  Section  2.5) 


Zi0.  +  LA*  +  =  0, 

(149) 

0*  +  y  =  0, 

(150) 

u,  +  +  y)  =  0, 

(151) 

+  G(0,-,0„u>i,  M  =  0. 

(152) 

The  vectors  of  streamfunction  and  vorticity  are  ordered  in  the  following  way 


where  the  i  subscript  denotes  evaluation  at  internal  nodes  while  the  s  subscript  de¬ 
notes  evaluation  at  the  airfoil  surface.  This  ordering  is  performed  to  isolate  equations 
that  explicitly  contain  time  derivatives.  The  linear  operators  Zj  and  L2  represent 
discrete  forms  of  the  Laplacian  operator  in  the  Poisson  equation  for  the  streamfunc¬ 
tion;  L3  is  the  discrete  form  of  the  operator  applied  at  the  airfoil  surface  for  the 
vorticity. 

The  temporal  development  of  small  perturbations  about  a  known,  steady-  state 
solution  (0°,w°)  is  assessed  by  making  the  collective  ansatz 

Mt)  =  0,°  +  e/<eAl,  Mi)  =  0.°  +  c/.cAt,  (153) 

Mt)  =  w®  +  e0,eA\  M1)  =  v°  +  tg,eXt.  (154) 

Substitution  into  Eq.  (150)  establishes  that  /,  =  0  and  further  substitution  into 
Eqs.  (149)  and  (151)  provides 


Lxii  +  9<  =  0, 


(155) 
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0. 


(156) 


L^fi  +  g,  = 

The  dynamic  vorticity  equation  yields  an  equation,  through  linearization,  for 
the  growth  coefficient  A 


G*/i  +  (<%  +  A)$,  +  Gl.9.  =  0,  (157) 

where  the  superscript  “o”  denotes  evaluation  at  (0°,o;0).  By  expressing  <7,  and  <7,  in 
terms  of  /,,  a  single  equation  for  /,  can  be  derived 


G%Ji  +  ( G cm  +  A)(-I1/,)  +  Gt.i-Lzfi)  =  0.  (158) 


This  equation  may  be  recast  as 


L?Afi  =  A/„  (159) 

where 

4  =  (GJ,  -  g-mlx  -  Gl'L,).  (160) 

The  method  of  Griewank  and  Reddien  may  now  be  applied  to  the  complete 
system  of  equations,  including  the  Hopf  point  criteria.  The  bifurcation  parameter, 
/?,  is  the  inverse  of  Reynolds  number  and  now  represents  an  unknown  variable.  The 
complete  system  is  then 


F\  —  +  Z'sV'*  +  <*>,•  =  0, 

(161) 

Fj  =  ij>.  +  y  =  0, 

(162) 

F3  =  L3(V>»  +  y)  +  w,  =  0, 

(163) 

F4  =  =  0, 

(164) 
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Fs  =  A(rl>i,xl>„u>i,w.;P)pi  +  6L1P2  =  0, 

(165) 

F6  =  A(rJ>i,  $)&  -  6L1P1  =  0, 

(166) 

F7  ~  qT Pi  =  0, 

(167) 

Fs  =  qTPi  -1  =  0, 

(168) 

where  p  =  pi  +  ipj  is  an  eigenvector  and  A  =  id  an  eigenvalue  satisfying  the  Hopf 
condition  for  Eq.  (159).  Eqs.  (164)  and  (166)  result  after  substituting  iO  into  Eq. 
(159)  and  equating  the  real  and  imaginary  parts.  Eqs.  (167)  and  (168)  represent 
two  real  normalizing  conditions  for  the  vector  p,  where  q  is  an  arbitrary  vector  with 
the  same  dimension  as  p.  Solution  of  the  complete  system  is  accomplished  through 
the  application  of  Newton’s  method  and  Gaussian  elimination.  The  implementation 
of  Newton’s  method  requires  a  linearization  of  Eqs.  (161-168)  in  order  to  apply  the 
iterative  scheme 

Av A"x  =  Av{xv+l  -  x ")  =  -F(xv)  =  —Fv,  (169) 

where  Av  is  the  Jacobian  matrix,  evaluated  at  the  v-th  Newton  iterate: 

(170) 

Xv 

The  iteration  procedure  is  continued  until  Fv,  or  alternately  A^x,  becomes  suffi¬ 
ciently  small.  Solutions  have  been  obtained  for  the  case  of  the  stationary  airfoil  and 
have  been  verified  by  time  integration  and  examination  of  the  eigenvalue  spectrum 
of  L^1  A  near  the  Hopf  point  by  Beran  and  Lutton  [9], 

C.2  Hopf  Bifurcation  Algorithm:  Moving  Airfoil 

The  following  approach  allows  a  similar  calculation  for  the  case  of  a  rotating, 
translating  airfoil.  To  accommodate  the  motion  of  the  airfoil,  an  extended  ansatz 
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set  is  required.  A  consistent  ansatz  set  would  begin,  in  the  case  of  the  pitch  axis, 
with  a  perturbation  about  the  angle  of  attack,  a  =  a0  -f  efaeXt.  However,  there 
are  several  disadvantages  in  such  an  approach.  Firstly,  the  surface  quantity  /,  is 
no  longer  zero,  and  the  complexity  and  memory  requirements  of  the  algorithm  are 
significantly  increased.  In  addition,  the  resulting  equation  set  is  no  longer  cast  in 
the  form  of  a  standard  eigenvalue  problem,  as  in  Eq.  (159),  invalidating  the  assumed 
approach.  A  brief  development  is  presented  to  indicate  the  difficulties.  The  following 
ansatz  set  is  chosen 

h  =  h°  +  efhex\  h  =  h°  +  eXfhext,  h  =  h°  +  eX2fheXt,  (171) 

a  =  a°  +  efaext,  a  =  a°  +  eXfaext,  a  =  a0  +  eX2faeXt.  (172) 

Application  of  Eqs.  (171)-(172)  to  Eqs.  (149)-(152)  provides  for  the  relationships 

Lift  +  Lif»  +  Q%  =  0,  (173) 

/.  +  yaf*  =  0,  (174) 

Lz(fi  +  yfa)  +  g$  =  0,  (175) 

+  G%Jt  +  G°w,g.  +  (A  +  G°Ui)9i  +  G° \fa  +  G°&\2fa  +  G-hX2fh  =  0.  (176) 

A  further  relationship  is  established  between  /<,/*,  and  fa  by  considering  the  airfoil 
equations  of  motion,  with  no  damping  or  structural  coupling  (cf.  Eqs.  (18)-(19)). 
Denoting  Qh  as  R  and  Qa  as  Q,  the  equation  representing  the  vertical  axis  may  be 
written  as 

mh  +  k^h  =  R. 

Linearization  provides 

(mAJ  +  h)U  =  (flj,  -  KM  -  KM)U 
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(177) 


-(K.y*  +  KMy*  +  K.L\y°)f<>, 

Similarly,  for  the  pitch  axis 

la  +  kaa  =  g, 

(/A3  +  k.  +  <%.$„  +  <&>$„  +  <&,£»„)/„  =  (Q*.  -  <£.£3  -  <3°,L,)/j.  (178) 

Returning  to  the  vorticity-transport  equation  and  substituting  provides 

-  A[£,  +  +  K  +  <?*.».  +  +  <&,£».)-'(<&,  -  Ql.U  -  (£.£,)]/, 

-  gj.(/a3  +  *„  +  Q;.y„  +  Ql.hy-  4  -  Ql, M  ~ 

-  czm-v a1 + *. + <?*.»« + + Q:M«r'(Q°*.  -  q:m  - 

-  G°UiL,ya(l*  +  a.  +  Q*.y»  +  <K.£3»„  +  <K,I3y„)-'(Q^  -  %  £3  -  Ql.L,)f. 

+  G“A(/A3  +  *„  +  Ql_y„  +  <£.£,»„  +  <£,£,»„)-■(%  ~  QtM  ~  )f- 

+  G|A3(/A3  +  *„  +  %$„  +  Ql'Lsy.  +  %  £»„)-'(<?;  -  Ql.U  -  <££,)/) 

+  GJA3(n>A3  +  *»)''[(«;,  -  «£.£3  -  <£■)  -  (**.»<■  +  7C.i3iG  +  /C.£.v») 

(/aj + 1„ + -  qi.l 3  -  <?:,£.)]/, 

+  (Gl,-Gl.L,-Gl:Ll)f.  =  0.  (179) 

This  equation  corresponds  to  Eq.  (159),  arrived  at  in  the  case  of  the  fixed  airfoil. 
However,  the  form  presented  does  not  allow  a  tractable  application  of  the  algorithm. 

To  rectify  these  problems,  a  hybrid  ansatz  set  is  applied: 

h  =  h°  +  efhext,  h  =  h°  +  t\fhtXi, 
a  =  a0  +  efhtXt,  a  =  a°  +  cA  fheXt.  (180) 

The  justification  for  this  approach  arises  from  the  examination  of  the  additional 
terms  appended  to  the  momentum  equation  (cf.  Eqs.  (9)  and  (10)).  The  additive 

terms  f\  and  fi  do  not  contain  a  or  h,  but  rather  refer  only  to  the  linear  and  angular 
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velocities  and  accelerations,  hence  only  these  terms  are  perturbed.  The  failure  to 
incorporate  a  perturbation  about  a  and  k  otherwise  has  the  effect  of  lagging  the 
boundary  conditions. 

Proceeding  with  this  approach,  the  substantial  change  from  the  original  equa¬ 
tion  set,  Eqs.  (149-152),  occurs  in  the  vorticity-transport  equation,  where 


Wt  +  G(^j,t/)„u;i,u>4,Q!,Q!,M  =  0. 


(181) 


Expanding  G  about  the  equilibrium  state  produces 


G(« +  «/.«".•  ••)-<?(«.—)  +  ef.e»  +  ~  cf,eu  +  ~  eg.eu 

Oxf)g  o  OXPi  o  OUg  o 

9G  ,,  8G  ,  ,,  dG  .  ,  ,, 
+  <9'e  +3a  ‘{°C  +9i 

•  O  0  o 


+  ^  a/.e*1  +  0(t!). 
Oh.  „ 


(182) 


Upon  subtraction  of  the  equilibrium  solution  from  the  expanded  vorticity-transport 
equqtion  and  division  by  ceAt,  an  equation  corresponding  to  Eq.  (158)  is  obtained: 


G%J.  +  G%Ji  +  (G°Wt  +  X)9i  +  G^g,  +  (G°6  +  A  G%)f0  +  A  G\fa  =  0.  (183) 


As  before,  f»  =  0,  and  the  relationship  between  fi,  g,,  and  g,  is  provided  by  Eqs. 
(155)  and  (156).  Then  a  relationship  between  /*,  /*,,  and  fa  is  developed  is  accom¬ 
plished  by  again  considering  the  airfoil  equations  of  motion,  Eqs.  (18)  and  (19),  with 
damping  and  structural  coupling  omitted.  The  approach  is  to  first  consider  the  un¬ 
steady  equation  to  establish  the  relationship  for  q,  and  then  consider  perturbations 
about  the  equilibrium  equation  to  establish  the  relationship  for  a. 


(Am,  +  kh)  fk  -  (/$,  -KM-  K.  I* )/.  =  0, 

(A I*  +  ka)fa  -  (q;(  -  qim  -  QIM)U  =  o- 
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(184) 

(185) 


Upon  substituting  Eqs.  (184)-(185)  into  Eq.  (183),  the  identical  form  is  obtained 
for  the  eigenvalue  problem  (cf.  Eq.  (159)),  provided  that  the  matrix  A  is  redefined 
as 


A  =  (G^-Gl.Lt-CC'L  3) 

+  (£g|  +  £<?*)  (OJ.  -  Qi.Lt  -  QI,U) 

+  (^Gt  +  lhGl)  (K.  -  K.L,  -  KM)  ■  (186) 

The  derivatives  of  G  with  respect  to  h,  h,  a,  and  a  can  be  evaluated  explicitly.  The 
vorticity-transport  equation  is  (cf.  Eq.  (29)) 

G  =  vwt  +  tkjv-  ^V2u>  -  {fx  -  f2)  =  0,  (187) 


where 


fx  =  -  h  sin  a  +  2 va  +  ya-  xa2, 
f2  =  h  cos  a  -  2nd  —  xa  —  yd 2, 


fi  =  Ui(  +  Vvfivi 
h  =  Zxfu  +  Vxfirr 


The  additional  contribution  to  the  vorticity  transport  equation  is  expressed  as 

(/1  -  ft)  =  iv  -  xJn  +  ydn  -  yvfn)-  (iss) 


The  required  derivatives  G%,  G|,  G?,  and  G°  Me  then 


= 


dG 

da 
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(189) 


d  ■ 

“  -5S«-« 

=  -^-1  faMi*  -  xnhi  +  ft/2*  -  y*/*) 
=  -2i-1(x{t>,,  -  x„t>€  +  y„u{  -  y*u„), 


similarly 


0G 

da 
d 

=  -«<*-*> 

-  d 

=  -^-1  Qz&th*  ~  xnfn  +  y«/2*  -  y.*/^) 
=  -  xvy(  +  yvx(  -  y(xv) 

=  -J~X{2J) 

=  -2, 


(190) 


and,  for  the  vertical  axis, 


<n 


dG 

dh 


dh(fi  h) 


=  -f  xirrOft/i* -  xvfu  +  Vihv  -  vM 
cth 

=  — «/-1(x((— sina),,  -  xv(-  sina)(  +  y^cosa),,  -  y„(cosa)f) 
=  0. 


(191) 


Since  h  does  not  appear  in  f\  or  /2, 


q  =  o. 


(192) 
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The  application  of  the  modified  algorithm,  as  described,  produced  an  insignif¬ 
icant  change  in  the  predicted  location  of  the  Hopf  point.  The  critical  Reynolds 
number,  Re^n,  obtained  using  the  original  algorithm  (without  the  extended  ansatz 
set)  for  Grid  1  was  550.  The  modified  algorithm  produced  a  nearly  identical  result, 
R^crit  =  552.  The  only  pronounced  effect  that  the  modifications  produced  was  a 
much  smaller  radius  of  convergence.  The  initial  values  of  Ia  and  ka  had  to  be  made 
very  large  ((3(107)),  i.e.,  approximating  a  fixed  airfoil.  After  obtaining  a  converged 
solution,  the  values  of  each  parameter  were  reduced  by  approximately  one  order  of 
magnitude,  and  a  new  solution  obtained.  This  procedure  was  applied  until  a  solution 
with  IQ  =  ka  =  100  was  obtained.  Throughout  this  process,  the  change  in  the  Hopf 
point  remained  insignificant.  This  was  a  rather  disappointing  result,  the  basis  of 
which  lies  in  the  nature  of  the  algorithm.  When  applying  Newton’s  method,  it  is  re¬ 
quired  to  evaluate  the  resulting  system  of  equations  at  equilibrium.  At  equilibrium, 
there  is  no  difference  between  the  fixed  airfoil  and  one  which  can  potentially  move, 
hence  no  difference  is  observed  in  the  location  of  the  Hopf  point.  Additionally,  it  is 
hypothesized  that  the  structural  components  principally  effect  the  stable,  limit  cycle 
solutions  in  the  bistable  region,  therefore,  a  method  which  essentially  searches  the 
path  of  equilibrium  solutions  (as  per  Figure  13)  for  an  instability  cannot  capture 
these  solutions. 

It  should  be  added,  as  a  criticism,  that  the  modified  approach  can  justly  be  cri¬ 
tiqued  as  lacking  rigor,  particularly  in  the  assumptions  associated  with  the  extended 
ansatz  set.  However,  additional  evidence  of  a  relatively  straightforward  nature  can 
be  offered  to  support  the  arguments  presented.  Consider  the  time-dependent  co¬ 
ordinate  transformation  presented  in  Part  II  (Eq.  (69)).  The  modifications  to  the 
boundary  equations  and  vorticity-transport  equation  are  presented  there,  in  Eqs. 
(71)-(86).  Note,  however,  that  the  grid  speeds  xT  and  yT  and  the  metric  terms  (T 
and  T)r  are  all  zero  at  equilibrium.  An  examination  of  the  aforementioned  equations 
reveals  that  when  applied  to  the  equilibrium  system  (Eqs.  (155-162)),  all  unsteady 
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terms  vanish  and  the  equations  assume  a  form  identical  to  those  of  the  fixed  airfoil, 
and  as  in  the  previous  case,  the  only  potential  differences  arise  in  the  application  of 
the  ansatz  and  the  redefinition  of  the  A  matrix. 
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Appendix  D.  Computer  Codes  and  Resources 


The  following  hardware  resources  were  utilized,  in  order  of  decreasing  usage, 
to  perform  the  computational  portions  of  this  dissertation 

•  AFIT  Kubota  (Titan  3000)  Cluster 

•  AFIT  Convex  C  220 

•  AFIT  Stardent  ST-2000 

•  Cray  XMP  (Wright-Patterson) 

•  Cray2  (Kirtland) 

D.l  Software  Documentation 

This  document  was  prepared  using  BTgX  on  an  AFIT  Sim  SPARCstation  2.  A 
short  synopsis  of  the  computer  codes  employed  is  provided,  along  with  a  description 
of  the  subroutines  contained  therein.  Four  main  codes  were  employed: 

1.  U22  -  incompressible  code  (Part  I) 

2.  U23  -  incompressible  code  (Part  II) 

3.  U26  -  incompressible  code  (Part  II) 

4.  BWO  -  compressible  code  (Part  II) 

The  incompressible  code  U22  allows  motion  of  the  airfoil  with  two  degrees  of 
freedom,  pitch  and  vertical  motion,  and  was  used  for  Part  I  of  the  dissertation.  U23 
is  identical  to  U22  except  for  the  implementation  of  the  Baldwin- Lomax  turbulence 
model  and  the  application  of  »7-upwinding.  U26  allows  airfoil  motion  in  the  pitch 
axis  only,  applying  the  time-dependent  coordinate  transformation  described  in  Sec¬ 
tion  4.2.  This  version  was  used  for  the  correlation  with  Theodorsen’s  function  and 
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contains  the  Baldwin-Lomax  turbulence  model,  with  17-upwinding  eliminated.  BWO 
is  the  Beam  and  Wanning  code  authored  by  Visbal  [27],  described  in  Section  4.4. 

A  short  description  of  input /output  files  for  the  above  codes  and  a  short  syn¬ 
opsis  of  the  subroutines  employed  is  provided  below. 

U22/23/26 

•  Input  Files 

udriver  -  input  file  to  specify  ffe,  At,  etc. 

foilgrd  -  grid  file 

foilout  -  restart  file  for  rj)  and  u> 

foiledv  -  restart  file  for  eddy  viscosity  (U23/26  only) 

•  Output  Files 

foilout  -  restart  file  for  %l>  and  a ; 

foiledv  -  restart  file  for  eddy  viscosity  (U23/26  only) 

foilvor  -  vorticity 

foilstr  -  streamfunction 

avt  -  angle  of  attack  vs.  time 

clvt  -  lift  coefficient  vs.  time 

elk  -  lift  coefficient  (theoretical)  vs.  time  (U26  only) 
hvt  -  vertical  displacement  vs.  time 

•  Subroutines 

AEOM  -  driver  file  for  solving  airfoil  equations  of  motion 

BANSOL  -  linear  system  solver  for  banded  matrices 

DECOMP  -  decomposes  A  matrix  into  lower  and  upper  triangular  matrices 
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DERIVS  -  computes  functional  derivatives  for  Runge-Kutta  integration 

ENERGY/ERG2  -  computes  aerodynamic  coefficients 

G14CLC  -  calculation  of  G  operator 

LUSOLV  -  solves  LU  decomposed  system 

L12CLC  -  calculates  L\  and  L2  operators 

L3CLC  -  calculates  L3  operator 

METRIC  -  calculates  transformation  metrics 

N1CLC  -  calculates  N\  vector 

N2CLC  -  calculates  N2  vector 

N3CLC  -  calculates  N3  vector 

N4CLC  -  calculates  N4  vector 

PACKER  -  calculates  column  index  for  banded  matrices 
RKD/RK4  -  implements  Runge-Kutta  integration 
TURBU  -  implements  Baldwin- Lomax  turbulence  model 
VISCO  -  calculates  coefficients  of  Laplace  operator,  fa 

BWO 

•  Input  Files 

odata  -  input  file  to  specify  -Re,  At,  etc. 

ogrid  -  grid  file;  also  contains  field  data  when  performing  restart 

•  Output  Files 

outo  -  output  file  for  lift,  drag,  norms,  etc. 

solno  -  output  file  for  field  variables,  also  used  for  restarts  when  renamed  as 
‘ogrid’ 
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avt  -  angle  of  attack  vs.  time 

clvt  -  lift  coefficient  vs.  time 

elk  -  lift  coefficient  (theoretical)  vs.  time 

hvt  -  vertical  displacement  vs.  time 

•  Subroutines 

BNDRY  -  implements  boundary  conditions 

BTRIDX/Y  -  block  tridiagonal  solvers  for  x  and  y  sweeps 

CMATA  -  computes  Jacobian  matrix  A 

CMATB  -  computes  Jacobian  matrix  B 

CM  ATM  -  computes  Jacobian  matrix  M 

CMATN  -  computes  Jacobian  matrix  N 

CMXU  -  computes  dynamic  viscosity  using  Sutherland’s  formula 
DAMPEX  -  implements  fourth-order  explicit  damping  in  f  direction 
DAMPEY  -  implements  fourth-order  explicit  damping  in  rj  direction 
INITIA  -  initialization  routine,  computes  common  parameters 
LIFT/1  -  computes  aerodynamic  coefficients 
METRIC  -  calculates  transformation  metrics  and  Jacobian 
OUTPUT  -  writes  flow  field  data  to  output  file 
RHSV  -  computes  RHS  vector 

SPECR  -  computes  scaling  factor  for  spectral  damping  coefficient 
STEPX  -  performs  (  sweep 
STEPY  -  performs  Tj  sweep 

TMSTEP  -  computes  times tep  (when  applying  local  timesteping) 
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TURB  -  implements  Baldwin- Lomax  turbulence  model 


The  computer  codes,  and  the  associated  drivers,  restart  files,  and  grids,  axe 
archived  on  the  AFIT  Convex  220  in  the  following  directories 

U22  -  /home/tap/mlutton/59x25 
U23  -  /home/tmp/mlutton/59x25/TURB 
U26  -  /home/tap/mlutton/UMOD 
BWO  -  /home/tap/mlutton/BWO 


D.2  Algorithm  Performance 

A  comparison  of  computer  code  performance  is  shown  in  Table  x.  The  compu¬ 
tations  were  performed  on  the  Kubota  cluster.  The  codes  were  executed  using  vector 
optimization  only.  The  incompressible  code  U26  is  seen  to  be  more  efficient  for  course 
grids.  The  increase  in  the  bandwidth  of  the  linear  system  results  in  a  degradation 
of  efficiency  for  the  incompressible  codes  as  the  grid  size  increases,  and  thus  the 
compressible  code  is  more  efficient  for  finer  grids.  A  more  efficient  implementation 
of  the  incompressible  code  has  been  developed  by  Beran  [56]. 


Code 

Grid 

No.  nodes 

Timesteps 

CPUsec 

CPUsec/node/iter 

59x25 

1475 

4000 

872.6 

wffwtnmi 

BWO 

65x25 

1625 

4000 

1355.2 

WEBESSSM 

U26 

99x40 

3960 

2000 

1731.5 

rnmnnsm 

BWO 

125x50 

6250 

2000 

2610.3 

U26 

139x50 

6950 

2000 

3804.2 

HggjjgjjEZQHH 

BWO 

209x108 

22572 

2000 

8608.1 

■BSBEBHii 

U26 

179x60 

10740 

2000 

6024.8 

HZiSSQIififli 

Table  5.  Comparison  of  Computer  Code  Performance 
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